## geotools-gt2-users

 [Geotools-gt2-users] Q on Geotools projections From: Davis Ford - 2009-02-12 18:58:36 ```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 ```
 Re: [Geotools-gt2-users] Q on Geotools projections From: Jody Garnett - 2009-02-12 22:54:35 ```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 > Geotools-gt2-users@... > https://lists.sourceforge.net/lists/listinfo/geotools-gt2-users > ```
 Re: [Geotools-gt2-users] Q on Geotools projections From: Davis Ford - 2009-02-13 02:58:25 ```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 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 >> Geotools-gt2-users@... >> 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 ```
 Re: [Geotools-gt2-users] Q on Geotools projections From: Jody Garnett - 2009-02-13 03:45:00 Attachments: Message as HTML ```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 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 > 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 > >> Geotools-gt2-users@... > >> 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 > ```
 Re: [Geotools-gt2-users] Q on Geotools projections From: Davis Ford - 2009-02-13 14:38:57 ```Thanks Jody -- this is a big help. Regards, Davis On Thu, Feb 12, 2009 at 10:44 PM, Jody Garnett 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 > 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 >> 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 >> >> Geotools-gt2-users@... >> >> 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 ```
 [Geotools-gt2-users] Help, geometry exception From: Malm Paul - 2009-03-20 13:27:35 ```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 ```
 Re: [Geotools-gt2-users] Help, geometry exception From: Michael Bedward - 2009-03-21 10:25:19 ```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 ```