Q on Geotools projections

classic Classic list List threaded Threaded
7 messages Options
Reply | Threaded
Open this post in threaded view
|

Q on Geotools projections

Davis Ford
Hi, I have a simple point obtained from Google Earth and I want to
calculate a buffer around this in units of miles.  So, I think the
general pattern is like this:

Geometry source= new WKTReader().read("POINT(-82.90755596903085
42.40409951227155)");
CoordinateReferenceSystem sourceCRS = CRS.decode("EPSG:4326");
// what should I project it to?
CoordinateReferenceSystem targetCRS = CRS.decode("???????");
MathTransform transform = CRS.findMathTransform(sourceCRS, targetCRS, true);

Geometry target = JTS.transform(source, transform);

// buffer by 1 mile (meters)
Geometry buffer = target.buffer(1609.344);

// re-project
transform = CRS.findMathTransform(targetCRS, sourceCRS, true);
Geometry geometry = JTS.transform(buffer, transform);

I'm not exactly sure what CRS I should project it to.  If I know
exactly in the world where this is I could use that, but is there a
general EPSG that I can use where the units of distance are in meters,
so I can use the JTS buffer method correctly?

Also, I read on the Wiki that Google is not really WGS84 -- will this
cause issues or should I try to use the CRS posted on the Wiki to get
more accurate results:

http://docs.codehaus.org/display/GEOTDOC/08+Google+Maps+Projection

Thanks in advance,
Davis

------------------------------------------------------------------------------
_______________________________________________
Geotools-gt2-users mailing list
[hidden email]
https://lists.sourceforge.net/lists/listinfo/geotools-gt2-users
Reply | Threaded
Open this post in threaded view
|

Re: Q on Geotools projections

jody.garnett
I am going to try and quickly answer; but you will need others to fill
in the details:
- I often know a good equal area projection for the area I will be
working in (EPSG:3005 bc albers is a good CRS for working in western
canada for example).
- you could study how bc albers is defined and make yourself a
CoordinateReferenceSystem on the fly for the specific point you are
interested in

Once you have the point in an equal area projection you can use buffer
to get a polygon; and then transform this polygon back to the google
projection if you like (it will show up as an ellipse).

Jody

Davis Ford wrote:

> Hi, I have a simple point obtained from Google Earth and I want to
> calculate a buffer around this in units of miles.  So, I think the
> general pattern is like this:
>
> Geometry source= new WKTReader().read("POINT(-82.90755596903085
> 42.40409951227155)");
> CoordinateReferenceSystem sourceCRS = CRS.decode("EPSG:4326");
> // what should I project it to?
> CoordinateReferenceSystem targetCRS = CRS.decode("???????");
> MathTransform transform = CRS.findMathTransform(sourceCRS, targetCRS, true);
>
> Geometry target = JTS.transform(source, transform);
>
> // buffer by 1 mile (meters)
> Geometry buffer = target.buffer(1609.344);
>
> // re-project
> transform = CRS.findMathTransform(targetCRS, sourceCRS, true);
> Geometry geometry = JTS.transform(buffer, transform);
>
> I'm not exactly sure what CRS I should project it to.  If I know
> exactly in the world where this is I could use that, but is there a
> general EPSG that I can use where the units of distance are in meters,
> so I can use the JTS buffer method correctly?
>
> Also, I read on the Wiki that Google is not really WGS84 -- will this
> cause issues or should I try to use the CRS posted on the Wiki to get
> more accurate results:
>
> http://docs.codehaus.org/display/GEOTDOC/08+Google+Maps+Projection
>
> Thanks in advance,
> Davis
>
> ------------------------------------------------------------------------------
> _______________________________________________
> Geotools-gt2-users mailing list
> [hidden email]
> https://lists.sourceforge.net/lists/listinfo/geotools-gt2-users
>  


------------------------------------------------------------------------------
Open Source Business Conference (OSBC), March 24-25, 2009, San Francisco, CA
-OSBC tackles the biggest issue in open source: Open Sourcing the Enterprise
-Strategies to boost innovation and cut costs with open source participation
-Receive a $600 discount off the registration fee with the source code: SFAD
http://p.sf.net/sfu/XcvMzF8H
_______________________________________________
Geotools-gt2-users mailing list
[hidden email]
https://lists.sourceforge.net/lists/listinfo/geotools-gt2-users
Reply | Threaded
Open this post in threaded view
|

Re: Q on Geotools projections

Davis Ford
Thank you Jody -- I appreciate the answer.

I'd be interested in figuring out how to generate a CRS on the fly.  I
will investigate this more.  Any pointers / tips / references are
greatly appreciated.

Final quick question: is there a general projected CRS that would work
for North America?  I am still learning all the different EPSG.
http://spatialreference.org/ is great.  For instance, I am using
EPSG:2807 and it does what I want if I insert that in the code I
pasted, because my point is in Michigan, but I am wondering if there
is one that would work (albeit with less accuracy) for North America
in general.

Regards,
Davis

On Thu, Feb 12, 2009 at 5:54 PM, Jody Garnett <[hidden email]> wrote:

> I am going to try and quickly answer; but you will need others to fill in
> the details:
> - I often know a good equal area projection for the area I will be working
> in (EPSG:3005 bc albers is a good CRS for working in western canada for
> example).
> - you could study how bc albers is defined and make yourself a
> CoordinateReferenceSystem on the fly for the specific point you are
> interested in
>
> Once you have the point in an equal area projection you can use buffer to
> get a polygon; and then transform this polygon back to the google projection
> if you like (it will show up as an ellipse).
>
> Jody
>
> Davis Ford wrote:
>>
>> Hi, I have a simple point obtained from Google Earth and I want to
>> calculate a buffer around this in units of miles.  So, I think the
>> general pattern is like this:
>>
>> Geometry source= new WKTReader().read("POINT(-82.90755596903085
>> 42.40409951227155)");
>> CoordinateReferenceSystem sourceCRS = CRS.decode("EPSG:4326");
>> // what should I project it to?
>> CoordinateReferenceSystem targetCRS = CRS.decode("???????");
>> MathTransform transform = CRS.findMathTransform(sourceCRS, targetCRS,
>> true);
>>
>> Geometry target = JTS.transform(source, transform);
>>
>> // buffer by 1 mile (meters)
>> Geometry buffer = target.buffer(1609.344);
>>
>> // re-project
>> transform = CRS.findMathTransform(targetCRS, sourceCRS, true);
>> Geometry geometry = JTS.transform(buffer, transform);
>>
>> I'm not exactly sure what CRS I should project it to.  If I know
>> exactly in the world where this is I could use that, but is there a
>> general EPSG that I can use where the units of distance are in meters,
>> so I can use the JTS buffer method correctly?
>>
>> Also, I read on the Wiki that Google is not really WGS84 -- will this
>> cause issues or should I try to use the CRS posted on the Wiki to get
>> more accurate results:
>>
>> http://docs.codehaus.org/display/GEOTDOC/08+Google+Maps+Projection
>>
>> Thanks in advance,
>> Davis
>>
>>
>> ------------------------------------------------------------------------------
>> _______________________________________________
>> Geotools-gt2-users mailing list
>> [hidden email]
>> https://lists.sourceforge.net/lists/listinfo/geotools-gt2-users
>>
>
>



--
Zeno Consulting, Inc.
home: http://www.zenoconsulting.biz
blog: http://zenoconsulting.wikidot.com
p: 248.894.4922
f: 313.884.2977

------------------------------------------------------------------------------
Open Source Business Conference (OSBC), March 24-25, 2009, San Francisco, CA
-OSBC tackles the biggest issue in open source: Open Sourcing the Enterprise
-Strategies to boost innovation and cut costs with open source participation
-Receive a $600 discount off the registration fee with the source code: SFAD
http://p.sf.net/sfu/XcvMzF8H
_______________________________________________
Geotools-gt2-users mailing list
[hidden email]
https://lists.sourceforge.net/lists/listinfo/geotools-gt2-users
Reply | Threaded
Open this post in threaded view
|

Re: Q on Geotools projections

jody.garnett
I did it once when implementing an "auto" projection for WMS support (see org.geotools.referencing.crs.AUTOCRSAuthorityFactory.AUTO42004) . But the real thing is here rather than looking up a formal definition using CRS.decode("EPSG:3005") ... you would be wither generating some WKT and then parsing it; or making use of the referencing module factories to construct the object you want an element at a time.

The user guide has some examples:
- http://docs.codehaus.org/display/GEOTDOC/03+CoordinateReferenceSystem

But as to the specifics of what you want to generate; lets look at the WKT definition of an equal area projection together and see what can be figured out - here is the definition of 3005... grabbed out of uDig 1.1.1 (which uses the older epsg-properties module).

PROJCS["NAD83 / BC Albers",
  GEOGCS["NAD83",
    DATUM["North_American_Datum_1983",
      SPHEROID["GRS 1980", 6378137.0, 298.257222101, AUTHORITY["EPSG","7019"]],
      TOWGS84[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
      AUTHORITY["EPSG","6269"]],
    PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]],
    UNIT["degree", 0.017453292519943295],
    AXIS["Lon", EAST],
    AXIS["Lat", NORTH],
    AUTHORITY["EPSG","4269"]],
  PROJECTION["Albers_Conic_Equal_Area"],
  PARAMETER["central_meridian", -126.0],
  PARAMETER["latitude_of_origin", 45.0],
  PARAMETER["standard_parallel_1", 50.0],
  PARAMETER["false_easting", 1000000.0],
  PARAMETER["false_northing", 0.0],
  PARAMETER["standard_parallel_2", 58.5],
  UNIT["m", 1.0],
  AXIS["x", EAST],
  AXIS["y", NORTH],
  AUTHORITY["EPSG","3005"]]

So you can see the axis are in meters which is good; the projection is an albers conic equal area which is what we want; I would start fiddling with the central meridian and latitude of origion to line this thing up with your data. I would consider using udig; or your own view to visualize what is going on as you change the values; and once you understand what they do you can write some java code to reproduce your result.

Jody


On Fri, Feb 13, 2009 at 1:58 PM, Davis Ford <[hidden email]> wrote:
Thank you Jody -- I appreciate the answer.

I'd be interested in figuring out how to generate a CRS on the fly.  I
will investigate this more.  Any pointers / tips / references are
greatly appreciated.

Final quick question: is there a general projected CRS that would work
for North America?  I am still learning all the different EPSG.
http://spatialreference.org/ is great.  For instance, I am using
EPSG:2807 and it does what I want if I insert that in the code I
pasted, because my point is in Michigan, but I am wondering if there
is one that would work (albeit with less accuracy) for North America
in general.

Regards,
Davis

On Thu, Feb 12, 2009 at 5:54 PM, Jody Garnett <[hidden email]> wrote:
> I am going to try and quickly answer; but you will need others to fill in
> the details:
> - I often know a good equal area projection for the area I will be working
> in (EPSG:3005 bc albers is a good CRS for working in western canada for
> example).
> - you could study how bc albers is defined and make yourself a
> CoordinateReferenceSystem on the fly for the specific point you are
> interested in
>
> Once you have the point in an equal area projection you can use buffer to
> get a polygon; and then transform this polygon back to the google projection
> if you like (it will show up as an ellipse).
>
> Jody
>
> Davis Ford wrote:
>>
>> Hi, I have a simple point obtained from Google Earth and I want to
>> calculate a buffer around this in units of miles.  So, I think the
>> general pattern is like this:
>>
>> Geometry source= new WKTReader().read("POINT(-82.90755596903085
>> 42.40409951227155)");
>> CoordinateReferenceSystem sourceCRS = CRS.decode("EPSG:4326");
>> // what should I project it to?
>> CoordinateReferenceSystem targetCRS = CRS.decode("???????");
>> MathTransform transform = CRS.findMathTransform(sourceCRS, targetCRS,
>> true);
>>
>> Geometry target = JTS.transform(source, transform);
>>
>> // buffer by 1 mile (meters)
>> Geometry buffer = target.buffer(1609.344);
>>
>> // re-project
>> transform = CRS.findMathTransform(targetCRS, sourceCRS, true);
>> Geometry geometry = JTS.transform(buffer, transform);
>>
>> I'm not exactly sure what CRS I should project it to.  If I know
>> exactly in the world where this is I could use that, but is there a
>> general EPSG that I can use where the units of distance are in meters,
>> so I can use the JTS buffer method correctly?
>>
>> Also, I read on the Wiki that Google is not really WGS84 -- will this
>> cause issues or should I try to use the CRS posted on the Wiki to get
>> more accurate results:
>>
>> http://docs.codehaus.org/display/GEOTDOC/08+Google+Maps+Projection
>>
>> Thanks in advance,
>> Davis
>>
>>
>> ------------------------------------------------------------------------------
>> _______________________________________________
>> Geotools-gt2-users mailing list
>> [hidden email]
>> https://lists.sourceforge.net/lists/listinfo/geotools-gt2-users
>>
>
>



--
Zeno Consulting, Inc.
home: http://www.zenoconsulting.biz
blog: http://zenoconsulting.wikidot.com
p: 248.894.4922
f: 313.884.2977


------------------------------------------------------------------------------
Open Source Business Conference (OSBC), March 24-25, 2009, San Francisco, CA
-OSBC tackles the biggest issue in open source: Open Sourcing the Enterprise
-Strategies to boost innovation and cut costs with open source participation
-Receive a $600 discount off the registration fee with the source code: SFAD
http://p.sf.net/sfu/XcvMzF8H
_______________________________________________
Geotools-gt2-users mailing list
[hidden email]
https://lists.sourceforge.net/lists/listinfo/geotools-gt2-users
Reply | Threaded
Open this post in threaded view
|

Re: Q on Geotools projections

Davis Ford
Thanks Jody -- this is a big help.

Regards,
Davis

On Thu, Feb 12, 2009 at 10:44 PM, Jody Garnett <[hidden email]> wrote:

> I did it once when implementing an "auto" projection for WMS support (see
> org.geotools.referencing.crs.AUTOCRSAuthorityFactory.AUTO42004) . But the
> real thing is here rather than looking up a formal definition using
> CRS.decode("EPSG:3005") ... you would be wither generating some WKT and then
> parsing it; or making use of the referencing module factories to construct
> the object you want an element at a time.
>
> The user guide has some examples:
> - http://docs.codehaus.org/display/GEOTDOC/03+CoordinateReferenceSystem
>
> But as to the specifics of what you want to generate; lets look at the WKT
> definition of an equal area projection together and see what can be figured
> out - here is the definition of 3005... grabbed out of uDig 1.1.1 (which
> uses the older epsg-properties module).
>
> PROJCS["NAD83 / BC Albers",
>   GEOGCS["NAD83",
>     DATUM["North_American_Datum_1983",
>       SPHEROID["GRS 1980", 6378137.0, 298.257222101,
> AUTHORITY["EPSG","7019"]],
>       TOWGS84[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
>       AUTHORITY["EPSG","6269"]],
>     PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]],
>     UNIT["degree", 0.017453292519943295],
>     AXIS["Lon", EAST],
>     AXIS["Lat", NORTH],
>     AUTHORITY["EPSG","4269"]],
>   PROJECTION["Albers_Conic_Equal_Area"],
>   PARAMETER["central_meridian", -126.0],
>   PARAMETER["latitude_of_origin", 45.0],
>   PARAMETER["standard_parallel_1", 50.0],
>   PARAMETER["false_easting", 1000000.0],
>   PARAMETER["false_northing", 0.0],
>   PARAMETER["standard_parallel_2", 58.5],
>   UNIT["m", 1.0],
>   AXIS["x", EAST],
>   AXIS["y", NORTH],
>   AUTHORITY["EPSG","3005"]]
>
> So you can see the axis are in meters which is good; the projection is an
> albers conic equal area which is what we want; I would start fiddling with
> the central meridian and latitude of origion to line this thing up with your
> data. I would consider using udig; or your own view to visualize what is
> going on as you change the values; and once you understand what they do you
> can write some java code to reproduce your result.
>
> Jody
>
>
> On Fri, Feb 13, 2009 at 1:58 PM, Davis Ford <[hidden email]>
> wrote:
>>
>> Thank you Jody -- I appreciate the answer.
>>
>> I'd be interested in figuring out how to generate a CRS on the fly.  I
>> will investigate this more.  Any pointers / tips / references are
>> greatly appreciated.
>>
>> Final quick question: is there a general projected CRS that would work
>> for North America?  I am still learning all the different EPSG.
>> http://spatialreference.org/ is great.  For instance, I am using
>> EPSG:2807 and it does what I want if I insert that in the code I
>> pasted, because my point is in Michigan, but I am wondering if there
>> is one that would work (albeit with less accuracy) for North America
>> in general.
>>
>> Regards,
>> Davis
>>
>> On Thu, Feb 12, 2009 at 5:54 PM, Jody Garnett <[hidden email]>
>> wrote:
>> > I am going to try and quickly answer; but you will need others to fill
>> > in
>> > the details:
>> > - I often know a good equal area projection for the area I will be
>> > working
>> > in (EPSG:3005 bc albers is a good CRS for working in western canada for
>> > example).
>> > - you could study how bc albers is defined and make yourself a
>> > CoordinateReferenceSystem on the fly for the specific point you are
>> > interested in
>> >
>> > Once you have the point in an equal area projection you can use buffer
>> > to
>> > get a polygon; and then transform this polygon back to the google
>> > projection
>> > if you like (it will show up as an ellipse).
>> >
>> > Jody
>> >
>> > Davis Ford wrote:
>> >>
>> >> Hi, I have a simple point obtained from Google Earth and I want to
>> >> calculate a buffer around this in units of miles.  So, I think the
>> >> general pattern is like this:
>> >>
>> >> Geometry source= new WKTReader().read("POINT(-82.90755596903085
>> >> 42.40409951227155)");
>> >> CoordinateReferenceSystem sourceCRS = CRS.decode("EPSG:4326");
>> >> // what should I project it to?
>> >> CoordinateReferenceSystem targetCRS = CRS.decode("???????");
>> >> MathTransform transform = CRS.findMathTransform(sourceCRS, targetCRS,
>> >> true);
>> >>
>> >> Geometry target = JTS.transform(source, transform);
>> >>
>> >> // buffer by 1 mile (meters)
>> >> Geometry buffer = target.buffer(1609.344);
>> >>
>> >> // re-project
>> >> transform = CRS.findMathTransform(targetCRS, sourceCRS, true);
>> >> Geometry geometry = JTS.transform(buffer, transform);
>> >>
>> >> I'm not exactly sure what CRS I should project it to.  If I know
>> >> exactly in the world where this is I could use that, but is there a
>> >> general EPSG that I can use where the units of distance are in meters,
>> >> so I can use the JTS buffer method correctly?
>> >>
>> >> Also, I read on the Wiki that Google is not really WGS84 -- will this
>> >> cause issues or should I try to use the CRS posted on the Wiki to get
>> >> more accurate results:
>> >>
>> >> http://docs.codehaus.org/display/GEOTDOC/08+Google+Maps+Projection
>> >>
>> >> Thanks in advance,
>> >> Davis
>> >>
>> >>
>> >>
>> >> ------------------------------------------------------------------------------
>> >> _______________________________________________
>> >> Geotools-gt2-users mailing list
>> >> [hidden email]
>> >> https://lists.sourceforge.net/lists/listinfo/geotools-gt2-users
>> >>
>> >
>> >
>>
>>
>>
>> --
>> Zeno Consulting, Inc.
>> home: http://www.zenoconsulting.biz
>> blog: http://zenoconsulting.wikidot.com
>> p: 248.894.4922
>> f: 313.884.2977
>
>



--
Zeno Consulting, Inc.
home: http://www.zenoconsulting.biz
blog: http://zenoconsulting.wikidot.com
p: 248.894.4922
f: 313.884.2977

------------------------------------------------------------------------------
Open Source Business Conference (OSBC), March 24-25, 2009, San Francisco, CA
-OSBC tackles the biggest issue in open source: Open Sourcing the Enterprise
-Strategies to boost innovation and cut costs with open source participation
-Receive a $600 discount off the registration fee with the source code: SFAD
http://p.sf.net/sfu/XcvMzF8H
_______________________________________________
Geotools-gt2-users mailing list
[hidden email]
https://lists.sourceforge.net/lists/listinfo/geotools-gt2-users
Reply | Threaded
Open this post in threaded view
|

Help, geometry exception

Paulm
Hi,
I'm using Geotools 2.5 for clipping maps into smaller tiles. I've downloaded a Vmap Level 0 map, eurasia.

when I shall clip the polygons with the following command: Geometry newGeom = featureGeometry.intersection(clipFeatureGeometry);
I'm getting the folloing error:
Exception in thread "main" com.vividsolutions.jts.geom.TopologyException: side location conflict [ (36.7401657104492, 45.4868316650391, NaN) ]
        at com.vividsolutions.jts.geomgraph.EdgeEndStar.propagateSideLabels(EdgeEndStar.java:298)
        at com.vividsolutions.jts.geomgraph.EdgeEndStar.computeLabelling(EdgeEndStar.java:136)
        at com.vividsolutions.jts.geomgraph.DirectedEdgeStar.computeLabelling(DirectedEdgeStar.java:127)
        at com.vividsolutions.jts.operation.overlay.OverlayOp.computeLabelling(OverlayOp.java:373)
        at com.vividsolutions.jts.operation.overlay.OverlayOp.computeOverlay(OverlayOp.java:173)
        at com.vividsolutions.jts.operation.overlay.OverlayOp.getResultGeometry(OverlayOp.java:127)
        at com.vividsolutions.jts.operation.overlay.OverlayOp.overlayOp(OverlayOp.java:66)
        at com.vividsolutions.jts.operation.overlay.snap.SnapOverlayOp.getResultGeometry(SnapOverlayOp.java:68)
        at com.vividsolutions.jts.operation.overlay.snap.SnapOverlayOp.overlayOp(SnapOverlayOp.java:25)
        at com.vividsolutions.jts.operation.overlay.snap.SnapIfNeededOverlayOp.getResultGeometry(SnapIfNeededOverlayOp.java:76)
        at com.vividsolutions.jts.operation.overlay.snap.SnapIfNeededOverlayOp.overlayOp(SnapIfNeededOverlayOp.java:25)
        at com.vividsolutions.jts.geom.Geometry.intersection(Geometry.java:1117)

When I look at the data in uDig I can see that the ocean area (Black sea) has a hole that isn't really a hole. E.g the hole starts in p1 goes to p2 and p3 then back to p2 and ends in p1. I.e the hole has no area.
This kind of holes I have also seen in DNC charts e.g a1316120 and others. I do not know why NGA creates those kinds of holes but it seems that it is common in VPF maps/DNC charts.

I have imported the VPF/DNC maps into a postgis database and clips the mapfeatures from those geometries. It is ok to import the data to the database and view them in uDig and QGIS, but it is not possible to cut them.

Does someone have a solution on this problem? This a real show stopper for me.


Kind regards,
Paul
------------------------------------------------------------------------------
Apps built with the Adobe(R) Flex(R) framework and Flex Builder(TM) are
powering Web 2.0 with engaging, cross-platform capabilities. Quickly and
easily build your RIAs with Flex Builder, the Eclipse(TM)based development
software that enables intelligent coding and step-through debugging.
Download the free 60 day trial. http://p.sf.net/sfu/www-adobe-com
_______________________________________________
Geotools-gt2-users mailing list
[hidden email]
https://lists.sourceforge.net/lists/listinfo/geotools-gt2-users
Reply | Threaded
Open this post in threaded view
|

Re: Help, geometry exception

mbedward
Hi Paul,

I don't have a solution (sorry) but there has been discussion of
dealing with very similar problems on the JTS list recently.  Have a
look at the list archives...

http://n2.nabble.com/jts-devel-f219725.html

...and/or pose your question to the folks there.

Hope this helps
Michael

------------------------------------------------------------------------------
Apps built with the Adobe(R) Flex(R) framework and Flex Builder(TM) are
powering Web 2.0 with engaging, cross-platform capabilities. Quickly and
easily build your RIAs with Flex Builder, the Eclipse(TM)based development
software that enables intelligent coding and step-through debugging.
Download the free 60 day trial. http://p.sf.net/sfu/www-adobe-com
_______________________________________________
Geotools-gt2-users mailing list
[hidden email]
https://lists.sourceforge.net/lists/listinfo/geotools-gt2-users