[
Date Prev][
Date Next][
Thread Prev][
Thread Next][
Date Index][
Thread Index]
[
List Home]
|
[udig-devel] Re: First (of many) reproject question
|
To Jody,
I did some tests and agree that this is weird. My test code follows (you
need to make a method in CRSService public to run it). I think the
problem is coming from CRSServices, since CTS is calculating the correct
points.
The following is my output. Notice how the second "Poly source" (coming
from poly1) is using NAD83 coordinates instead of BC Albers. This
indicates to me that CRSServices is somehow modifying its input
geometries. However, when transforming the envelopes, poly1 is back to
its original self (spooky).
Source CRS: PROJCS["BC_Albers", GEOGCS["GCS_North_American_1983",
DATUM["D_North_American_1983", SPHEROID["GRS_1980", 6378137.0,
298.257222101], TOWGS84[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]],
PRIMEM["Greenwich", 0.0], UNIT["degree of angle",0.017453292519943295],
AXIS["Longitude",EAST], AXIS["Latitude",NORTH]],
PROJECTION["Albers_Conic_Equal_Area"], PARAMETER["semi_major",
6378137.0], PARAMETER["semi_minor", 6356752.314140356],
PARAMETER["central_meridian", -126.0], PARAMETER["latitude_of_origin",
45.0], PARAMETER["false_easting", 1000000.0],
PARAMETER["false_northing", 0.0], PARAMETER["standard_parallel_1",
50.0], PARAMETER["standard_parallel_2", 58.5], UNIT["metre",1.0],
AXIS["x",EAST], AXIS["y",NORTH]]
Target CRS: GEOGCS["GCS_North_American_1983",
DATUM["D_North_American_1983", SPHEROID["GRS_1980", 6378137.0,
298.257222101], TOWGS84[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]],
PRIMEM["Greenwich", 0.0], UNIT["degree of angle",0.017453292519943295],
AXIS["Longitude",EAST], AXIS["Latitude",NORTH]]
Math Trans: INVERSE_MT[PARAM_MT["Albers_Conic_Equal_Area",
PARAMETER["semi_major",6378137.0],
PARAMETER["semi_minor",6356752.314140356],
PARAMETER["central_meridian",-126.0],
PARAMETER["latitude_of_origin",45.0],
PARAMETER["false_easting",1000000.0], PARAMETER["false_northing",0.0],
PARAMETER["standard_parallel_1",50.0],
PARAMETER["standard_parallel_2",58.5]]]
Poly source :POLYGON ((1187128 395268, 1187128 396027, 1188245 396027,
1188245 395268, 1187128 395268))
Env source :Env[1187128.0 : 1188245.0, 395268.0 : 396027.0]
Poly source :POLYGON ((-123.47009555832284 48.54326156200811,
-123.46972894676578 48.55009592211517, -123.45463828850829
48.54973520360885, -123.4550070827961 48.54290089163769,
-123.47009555832284 48.54326156200811))
Poly computed:POLYGON ((-123.47009555832284 48.54326156200811,
-123.46972894676578 48.55009592211517, -123.45463828850829
48.54973520360885, -123.4550070827961 48.54290089163769,
-123.47009555832284 48.54326156200811))
Poly expected:POLYGON ((-123.470095558323 48.5432615620081,
-123.469728946766 48.5500959221152, -123.454638288508 48.5497352036088,
-123.455007082796 48.5429008916377, -123.470095558323 48.5432615620081))
Env source :Env[1187128.0 : 1188245.0, 395268.0 : 396027.0]
Env source :Env[1187128.0 : 1188245.0, 395268.0 : 396027.0]
Env computed:Env[-123.4700927734375 : -123.45463562011719,
48.54290008544922 : 48.55009460449219]
Env expected:Env[-123.470095558323 : -123.454638288508,
48.5429008916377 : 48.5500959221152]
Now, when I comment out the three statements printing the first poly and
envelope, I get the following results. The source envelope that is being projected
is already in NAD83, giving the same wrong outputs you had in your test.
[WKT cut]
Poly source :POLYGON ((-123.47009555832284 48.54326156200811, -123.46972894676578 48.55009592211517, -123.45463828850829 48.54973520360885, -123.4550070827961 48.54290089163769, -123.47009555832284 48.54326156200811))
Poly computed:POLYGON ((-123.47009555832284 48.54326156200811, -123.46972894676578 48.55009592211517, -123.45463828850829 48.54973520360885, -123.4550070827961 48.54290089163769, -123.47009555832284 48.54326156200811))
Poly expected:POLYGON ((-123.470095558323 48.5432615620081, -123.469728946766 48.5500959221152, -123.454638288508 48.5497352036088, -123.455007082796 48.5429008916377, -123.470095558323 48.5432615620081))
Env source :Env[-123.47009555832284 : -123.45463828850829, 48.54290089163769 : 48.55009592211517]
Env source :Env[-123.47009555832284 : -123.45463828850829, 48.54290089163769 : 48.55009592211517]
Env computed:Env[-138.4474639892578 : -138.4474639892578, 44.199676513671875 : 44.199676513671875]
Env expected:Env[-123.470095558323 : -123.454638288508, 48.5429008916377 : 48.5500959221152]
Hopefully this is helpful. Good luck sorting this out.
Rueben
/*
* JodysCRSProblem.java
*
* Created on December 8, 2004, 8:38 PM
*/
import org.geotools.cs.CoordinateSystem;
import org.geotools.data.crs.CRSService;
import org.geotools.ct.MathTransform;
import com.vividsolutions.jts.geom.Polygon;
import com.vividsolutions.jts.geom.GeometryFactory;
import com.vividsolutions.jts.geom.LinearRing;
import com.vividsolutions.jts.geom.Coordinate;
import com.vividsolutions.jts.geom.Envelope;
/**
*
* @author rschulz
*/
public class JodysCRSProblem {
CRSService service = new CRSService();
GeometryFactory fact = new GeometryFactory();
/** Creates a new instance of JodysCRSProblem */
public JodysCRSProblem() throws Exception {
testTranform();
}
public void testTranform() throws Exception {
CoordinateSystem bc =
service.createCoordinateSystem("EPSG:42102");
CoordinateSystem latlong =
service.createCoordinateSystem("EPSG:4269");
MathTransform transform = CRSService.reproject( bc, latlong,
true );
System.out.println("Source CRS: " + bc);
System.out.println("Target CRS: " + latlong);
System.out.println("Math Trans: " + transform + "\n");
// origional bc alberts
Polygon poly1 = poly( new double[] {
1187128,395268, 1187128,396027,
1188245,396027, 1188245,395268,
1187128,395268} );
// transformed
Polygon poly2 = poly( new double[] {
-123.470095558323,48.5432615620081,
-123.469728946766,48.5500959221152,
-123.454638288508,48.5497352036088,
-123.455007082796,48.5429008916377,
-123.470095558323,48.5432615620081} );
//the following 3 statements were commented out in the second test.
System.out.println( "Poly source :"+ poly1 );
Envelope before = poly1.getEnvelopeInternal();
System.out.println("Env source :" + before + "\n");
Polygon polyAfter = CRSService.transform( poly1, transform );
System.out.println( "Poly source :"+ poly1 );
System.out.println( "Poly computed:"+ polyAfter );
System.out.println( "Poly expected:"+ poly2 + "\n");
//assertEquals( poly2, polyAfter );
before = poly1.getEnvelopeInternal();
System.out.println("Env source :" + before );
Envelope after = CRSService.transform( before, transform );
System.out.println("Env source :" + before );
System.out.println("Env computed:" + after );
System.out.println("Env expected:" + poly2.getEnvelopeInternal()
+ "\n");
//assertEquals( poly2.getEnvelopeInternal(), after );
}
Polygon poly (double[] pts) {
Coordinate[] points = new Coordinate[pts.length/2];
for(int i=0, j=0; i<pts.length; i=i+2, j++){
points[j] = new Coordinate(pts[i], pts[i+1], 0.0);
}
LinearRing ring = fact.createLinearRing(points);
return fact.createPolygon(ring, null);
}
/**
* @param args the command line arguments
*/
public static void main(String[] args) throws Exception {
// TODO code application logic here
new JodysCRSProblem();
}
}
On Wed, 2004-08-12 at 11:46 -0800, Jody Garnett wrote:
> Hi CRS guys I have a question / testcase for you ... the code is at the
> end of this message.
>
> In short I am getting very different results from the Postgis
> transformation code - I am sure we went over this
> back in the spring. I recall at the end of the day we got both systems
> returning with in 10cm of each other so
> I am probably doing something wrong. The epsg.properties file being used
> was generated from PostGIS so
> the WKT definitions they are working from should be identical.
>
> Now the answer may be one of:
> 0) I messed up
> 1) Switch over to the CRS system
> 2) The old CS objects no longer work, or parse WKT differently
> 4) The epsg.properties file is some how wrong
>
> If I messed up on this one I really would like to fix things as soon as
> possible because all DataStores, Lite renderer and UDIG
> are funneling through the CRSService API. It may be that CRSService
> wants to go and live next to all the other CRS code,
> it operates as a Facade for that module after all....
>
> Jody Garnett
>
> ------------------------
> public void testTranform() throws Exception {
> CoordinateSystem bc = service.createCoordinateSystem("EPSG:42102");
> CoordinateSystem latlong =
> service.createCoordinateSystem("EPSG:4269");
>
> MathTransform transform = CRSService.reproject( bc, latlong, true );
>
> // origional bc alberts
> Polygon poly1 = poly( new double[] {
> 1187128,395268, 1187128,396027,
> 1188245,396027, 1188245,395268,
> 1187128,395268} );
>
> // transformed
> Polygon poly2 = poly( new double[] {
> -123.470095558323,48.5432615620081,
> -123.469728946766,48.5500959221152,
> -123.454638288508,48.5497352036088,
> -123.455007082796,48.5429008916377,
> -123.470095558323,48.5432615620081} );
>
> Polygon polyAfter = CRSService.transform( poly1, transform );
> System.out.println( "expected:"+ polyAfter );
> System.out.println( "expected:"+ poly2 );
> //assertEquals( poly2, polyAfter );
>
> Envelope before = poly1.getEnvelopeInternal();
> Envelope after = CRSService.transform( before, transform );
> assertEquals( poly2.getEnvelopeInternal(), after );
> }
> ----------------------
> Here is the output:
> actual:POLYGON ((-123.47009555832284 48.54326156200811,
> -123.46972894676578 48.55009592211517, -123.45463828850829
> 48.54973520360885, -123.4550070827961 48.54290089163769,
> -123.47009555832284 48.54326156200811))
> expected:POLYGON ((-123.470095558323 48.5432615620081, -123.469728946766
> 48.5500959221152, -123.454638288508 48.5497352036088, -123.455007082796
> 48.5429008916377, -123.470095558323 48.5432615620081))
>
> actual:Env[-138.4474639892578 : -138.4474639892578, 44.199676513671875
> : 44.199676513671875]
> expected:Env[-123.470095558323 : -123.454638288508, 48.5429008916377 :
> 48.5500959221152]
