Andrea Aime ha scritto:
> Andrea Aime ha scritto:
>
>> Now, I'm wondering, besides extending point sampling inside
>> the original bbox (instead of using only the perimeter) do
>> we have any smarter method? Even sampling inside the bbox
>> we could loose part of the actual map area in 4326 if sampling
>> did not happen to catch the pole...
Ok, I rolled a solution, which is basically a brute force recursive
sampling, quadtree inspired, that ends its recursion when the area
of the current envelope did not expand more than a certain tolerance.
Seems to be quick enough, and to converge relatively fast to the
ideal envelope even in this pesky stereographic case.
Here is a code sample:
-------------------------------------------------------------
import java.util.Iterator;
import org.geotools.geometry.jts.JTS;
import org.geotools.geometry.jts.ReferencedEnvelope;
import org.geotools.referencing.CRS;
import org.opengis.geometry.Envelope;
import org.opengis.referencing.crs.CoordinateReferenceSystem;
import org.opengis.referencing.operation.MathTransform;
import org.opengis.referencing.operation.TransformException;
public class ReprojectEnvelopeTest {
public static void main(String[] args) throws Exception {
CoordinateReferenceSystem c3031 = CRS.decode("EPSG:3031", true);
CoordinateReferenceSystem c4326 = CRS.decode("EPSG:4326", true);
ReferencedEnvelope envelope = new
ReferencedEnvelope(-3943612.4042124213,
3729092.5890516187, -4078471.954436003,
4033483.085688618, c3031);
MathTransform transform = CRS.findMathTransform(c3031, c4326);
Envelope out = CRS.transform(transform, envelope);
System.out.println("CRS reproject: " + out);
envelope = new ReferencedEnvelope(-3943612.4042124213,
3729092.5890516187,
-4078471.954436003, 4033483.085688618, c3031);
System.out.println("JTS transform: " + JTS.transform(envelope,
transform));
System.out.println("Recursive sample");
envelope = new ReferencedEnvelope(-3943612.4042124213,
3729092.5890516187,
-4078471.954436003, 4033483.085688618, c3031);
double tolerance = 0.1;
for (int i = 0; i < 20; i++) {
transform(envelope, transform, tolerance);
tolerance /= 10.0;
}
}
public static com.vividsolutions.jts.geom.Envelope transform(final
ReferencedEnvelope envelope, final MathTransform transform, double
tolerance) throws TransformException
{
System.out.println("Tolerance: " + tolerance);
long start = System.currentTimeMillis();
com.vividsolutions.jts.geom.Envelope result = new
com.vividsolutions.jts.geom.Envelope();
// initial sampling, 4 corners of the envelope
final double xmin = envelope.getMinX();
final double xmax = envelope.getMaxX();
final double ymin = envelope.getMinY();
final double ymax = envelope.getMaxY();
double[] coords = new double[8];
coords[0] = xmin;
coords[1] = ymax;
coords[2] = xmax;
coords[3] = ymax;
coords[4] = xmin;
coords[5] = ymin;
coords[6] = xmax;
coords[7] = ymin;
transform.transform(coords, 0, coords, 0, 4);
for (int i = 0, j = 0; i < 4; i++) {
result.expandToInclude(coords[j++], coords[j++]);
}
// recursive quadtree scan, each call samples a cross shaped
set of 5 point
int samples = transform(envelope.getCenter(0),
envelope.getCenter(1), envelope.getWidth() / 2, envelope.getHeight() /
2, result, transform, new double[10], tolerance);
System.out.println("Elapsed " + (System.currentTimeMillis() -
start) + ", " + samples + " samples " + result );
return result;
}
private static int transform(double cx, double cy, double hw,
double hh, com.vividsolutions.jts.geom.Envelope result, MathTransform
transform, double[] coords, double tolerance) throws TransformException {
double initialArea = result.getWidth() * result.getHeight();
// build the cross, first center, then 4 satellites
coords[0] = cx;
coords[1] = cy;
coords[2] = cx;
coords[3] = cy - hh;
coords[4] = cx;
coords[5] = cy + hh;
coords[6] = cx - hw;
coords[7] = cy;
coords[8] = cx + hw;
coords[9] = cy;
transform.transform(coords, 0, coords, 0, 5);
for (int i = 0, j = 0; i < 4; i++) {
result.expandToInclude(coords[j++], coords[j++]);
}
// if we did expand the envelope at least one thousands, keep
on recursinv it
int samples = 10;
if(result.getWidth() * result.getHeight() / initialArea > (1 +
tolerance)) {
double qw = hw / 2;
double qh = hh / 2;
samples += transform(cx - qw, cy - qh, qw, qh, result,
transform, coords, tolerance);
samples +=transform(cx + qw, cy - qh, qw, qh, result,
transform, coords, tolerance);
samples +=transform(cx - qw, cy + qh, qw, qh, result,
transform, coords, tolerance);
samples +=transform(cx + qw, cy + qh, qw, qh, result,
transform, coords, tolerance);
}
return samples;
}
}
----------------------------------------------------------------------
If you run it, you'll see that with a tolerance of 0.0001 we're already
pretty much within the ideal envelope, and even asking for a tolerance
of 1e-20 does not loads the processor up enough to make for a measurable
delay (the sysouts printing elapsed time are more or less always
0 on my PC, with some occasional 15 milliseconds, I guess this is the
shortest time my PC tick can measure).
Cheers
Andrea
-------------------------------------------------------------------------
This SF.net email is sponsored by DB2 Express
Download DB2 Express C - the FREE version of DB2 express and take
control of your XML. No limits. Just data. Click to get it now.
http://sourceforge.net/powerbar/db2/
_______________________________________________
Geotools-devel mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/geotools-devel
