tag:blogger.com,1999:blog-80427704415021490162024-02-22T06:31:17.172-08:00GeoScriptJustin Deoliveirahttp://www.blogger.com/profile/04340292357363191337noreply@blogger.comBlogger13125tag:blogger.com,1999:blog-8042770441502149016.post-22114101905059922072011-07-23T13:24:00.001-07:002011-07-24T21:50:18.142-07:00GeoScript at FOSS4G 2011The GeoScript project has been quiet for a while now, but that does not mean developers haven't been busy. With the <a href="http://2011.foss4g.org/">FOSS4G</a> conference approaching development is ramping up. GeoScript will have a good showing this year, with <a href="http://2011.foss4g.org/sessions/geoscript-spatial-capabilities-scripting-languages">two talks</a> and one <a href="http://2011.foss4g.org/sessions/tutorial-working-geoscript">tutorial</a> on the conference program.<div><br /></div><div>The first talk, being delivered by Jared Erickson and myself, is titled <a href="http://2011.foss4g.org/sessions/geoscript-spatial-capabilities-scripting-languages">GeoScript - Spatial Capabilities for Scripting Languages</a> and will be an introduction of the project as a whole. The second talk titled <a href="http://2011.foss4g.org/sessions/scripting-geoserver-geoscript">Scripting GeoServer with GeoScript</a> will be presented by Tim Schaub and myself and will focus in more on GeoScript as a scripting engine integrated with GeoServer.</div><div><br /></div><div>The tutorial <a href="http://2011.foss4g.org/sessions/tutorial-working-geoscript">Working with GeoScript</a> will be a 1.5 hour session and will provide a working introduction to GeoScript that will allow attendees to get their hands dirty with some code. </div><div><br /></div><div>Developer turn out is going to be great this year. All four of David, Jared, Tim and myself will be attending this year so the four primary implementations will all be well represented. As well Ivan Willig is currently working hard on an implementation of <a href="https://github.com/iwillig/geoscript-clj">GeoScript for Clojure</a> that looks very interesting. Ivan hopes to get something up and on the site by FOSS4G.</div><div><br /></div><div>Can't wait for the conference and look forward to seeing everyone there!</div>Justin Deoliveirahttp://www.blogger.com/profile/04340292357363191337noreply@blogger.com1tag:blogger.com,1999:blog-8042770441502149016.post-8406990175239302312011-03-04T18:30:00.001-08:002011-08-02T07:55:20.142-07:00Styling with GeoScriptOver the past few months one of the major developments in GeoScript has been an API for styling. I finally got around to implementing the api for GeoScript python and am happy with how things are progressing.<br /><br />GeoScript, being based on GeoTools, naturally uses SLD as the underlying styling engine. Those who have used SLD know that it is a complicated beast. While it is a very powerful symbology language, it is also very complex and verbose. And for those not familiar with SLD and XML it represents a substantial learning curve. GeoScript hopes to make SLD a bit more approachable.<br /><br />Let's start off with a simple example:<br /><br /><pre class="brush:py">from geoscript.geom import *<br />from geoscript.style import *<br />from geoscript.render import draw<br /><br />style = Stroke('black', 2) + Fill('gray', 0.75)<br />g = Point(0,0).buffer(0.2)<br /><br />draw(g, style, (250,250))<br /></pre><br /><br />Which results in:<br /><br /><img style="width: 250px; height: 250px;" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEih6luVBZGDbiF1HTkWsrzrHJlgOulABkvaZwV70FD0L8nBx32idBEdADUxXT4VhOWfGcX0lrR2e3MLaEtSWp3A69qyeHwDK12-iZRGG_h2sEPgAn8SwIxVnVFwuW5j000Icwb1qibAxS66/s320/simple.png" border="0" alt="" id="BLOGGER_PHOTO_ID_5632951921944497314" /><br /><br />The equivalent SLD for this style is:<br /><pre class="brush:xml"><sld:UserStyle xmlns="http://www.opengis.net/sld" xmlns:sld="http://www.opengis.net/sld" xmlns:ogc="http://www.opengis.net/ogc"><br /><sld:FeatureTypeStyle><br /><sld:Name>My Style</sld:Name><br /><sld:Rule><br /><sld:LineSymbolizer><br /><sld:Stroke><br /><sld:CssParameter name="stroke-width">2</sld:CssParameter><br /></sld:Stroke><br /></sld:LineSymbolizer><br /><sld:PolygonSymbolizer><br /><sld:Fill><br /><sld:CssParameter name="fill-opacity">0.75</sld:CssParameter><br /></sld:Fill><br /></sld:PolygonSymbolizer><br /></sld:Rule><br /></sld:FeatureTypeStyle><br /></sld:UserStyle><br /></pre><br /><br />Which is not all that bad as far as SLD goes. Let's look at a more complex example, one that anyone who as ever used GeoServer will be familiar with. The all too well known <a href="http://demo.opengeo.org/geoserver/wms?service=WMS&version=1.1.0&request=GetMap&layers=topp:states&styles=&bbox=-124.731,24.956,-66.97,49.372&width=780&height=330&srs=EPSG:4326&format=application/openlayers">states</a> layer that is styled with the <a href="http://demo.opengeo.org/geoserver/styles/popshade.sld">popshade</a> SLD. As far as SLD styles go popshade is actually relatively simple, but still involves hacking out about 100 lines of XML. What does the equivalent style look like with GeoScript?<br /><br /><pre class="brush:py">style = Stroke()<br />style += Fill('#4DFF4D', 0.7).where('PERSONS < 2000000')<br />style += Fill('#FF4D4D', 0.7).where('PERSONS BETWEEN 2000000 AND 4000000')<br />style += Fill('#4D4DFF', 0.7).where('PERSONS > 4000000')<br />style += Label('STATE_ABBR').font('14 "Times New Roman"')<br /></pre><br /><br />A bit more concise than the SLD equivalent. But for those who just can't go without their SLD remember that this is all geotools styling under the hood, so getting at the SLD is easy. So:<br /><br /><pre class="brush:py">style.asSLD(sys.stdout)<br />style.asSLD('popshade.sld')<br /></pre><br /><br />Dumps out the SLD that looks very similar to popshade.sld above.<br /><br />Continuing on with another example of a style type that is frequently used: line casing. More specifically rendering two lines on top of each other with the "bottom" line styled a few pixels thicker than the "top" line to give the appearance of a border or casing. This technique is commonly used to style road data.<br /><br /><img style="width: 300px; height: 300px;" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEj49zT5TCMdg-568r_7qybIzD14Xf30Jyg8DPfSwfsMj1e5dNSgidd82Jkpf75iNySfqcY6vC59W5loJ5RfjB6i5teOKRtjVOqUB56WJemzmbP8rw3tVM-9v66Ly8kNaLO4imh9JmJWHKN0/s320/line.png" border="0" alt="" id="BLOGGER_PHOTO_ID_5634125753096721970" /><br /><br />First the SLD:<br /><br /><pre class="brush:xml"><FeatureTypeStyle><br /><Rule><br /><LineSymbolizer><br /><Stroke><br /><CSSParameter name="stroke">#000000</CSSParameter><br /><CSSParameter name="stroke-width">5</CSSParameter><br /></Stroke><br /></LineSymbolizer><br /></Rule><br /></FeatureTypeStyle><br /><FeatureTypeStyle><br /><Rule><br /><LineSymbolizer><br /><Stroke><br /><CSSParameter name="stroke">#ffffff</CSSParameter><br /><CSSParameter name="stroke-width">3</CSSParameter><br /></Stroke><br /></LineSymbolizer><br /></Rule><br /></FeatureTypeStyle><br /></pre><br /><br />Not bad, but in GeoScript it is a one liner:<br /><br /><pre class="brush:py">style = Stroke('#000000', 6) + Stroke('#ffffff', 3).zindex(1)<br /></pre><br /><br />Let's add some labels to those lines:<br /><br /><pre class="brush:py">style += Label('name', '12pt Arial').linear(follow=True,group=True,repeat=25).halo(Fill('#ffffff'), 2)<br /></pre><br /><br /><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg6NhMtH2i9y70-2hxv2dbZwEE9UklEhu5RA5IUV8TXMUsN26bEHDs8EB9qp8s5gvYYfsXj5YksCAXTKd4I1OQtaPCkwJzEfXnxA6q18JtRslI6hv4pXb6bAm-piZm56fdt0swoLQYk5kyN/s1600/line_label.png" onblur="try {parent.deselectBloggerImageGracefully();} catch(e) {}"><img style="cursor:pointer; cursor:hand;width: 320px; height: 320px;" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg6NhMtH2i9y70-2hxv2dbZwEE9UklEhu5RA5IUV8TXMUsN26bEHDs8EB9qp8s5gvYYfsXj5YksCAXTKd4I1OQtaPCkwJzEfXnxA6q18JtRslI6hv4pXb6bAm-piZm56fdt0swoLQYk5kyN/s320/line_label.png" border="0" alt="" id="BLOGGER_PHOTO_ID_5634127701299204690" /></a><br /><br />Those are just a few examples. Check out the <a href="http://geoscript.org/py/api/style/">api docs</a> for a complete list of what is available. More to come soon so stay tuned.Justin Deoliveirahttp://www.blogger.com/profile/04340292357363191337noreply@blogger.com0tag:blogger.com,1999:blog-8042770441502149016.post-76324652653515849322010-06-29T22:57:00.000-07:002010-07-03T17:09:49.369-07:00Unwrapped JTS with PythonWhile Justin has been <a href="http://github.com/jdeolive/geoscript-py">hard at work</a> wrapping <a href="http://www.vividsolutions.com/jts/jtshome.htm">JTS</a> and <a href="http://www.geotools.org/">GeoTools</a> in clear and simple Python syntax, I’ve been tinkering with some JTS classes that aren’t wrapped in GeoScript Python. Well, not wrapped yet. One of the great things about using Jython (or any scripting language written for the Java platform) is that you can still explore JTS and GeoTools even if some functionality isn’t implemented in GeoScript.<br /><br />Jared produced a convex hull and minimum bounding circle with GeoScript Groovy in a recent <a href="http://geoscriptblog.blogspot.com/2010/05/calculating-convex-hull-and-minimum.html">post</a>. I tried to replicate this workflow but quickly realized certain classes and functions weren’t available in GeoScript Python. Did it matter? Not very much in this case. Here’s the convex hull and minimum bounding circle for a shapefile of Washington weather stations produced in a previous <a href="http://geoscriptblog.blogspot.com/2010/06/saving-weather-stations-to-shapefile.html">post</a>. Pay attention to what’s being imported.<br /><pre class="brush:py"><br />from geoscript.layer import Shapefile<br />from geoscript.viewer import plot<br />from com.vividsolutions.jts.geom import GeometryCollection<br />from com.vividsolutions.jts.geom import GeometryFactory<br />from com.vividsolutions.jts.geom import PrecisionModel<br />from com.vividsolutions.jts.algorithm import MinimumBoundingCircle<br /><br />#Get the point shapefile<br />staLayer = Shapefile('/home/gregcorradini/wa_stations.shp')<br /><br />#Create a list of feature geometries<br />pntGeoms = [f.geom for f in staLayer.features()]<br /><br />#Create a GeometryCollection from our list of points<br />geomColl = GeometryCollection(pntGeoms,GeometryFactory(PreciscionModel(),4326))<br /><br />#Get the geometry collection's convex hull<br />geomConvexH = geomColl.convexHull()<br /><br />#Get the geometry collection's minimum bounding circle<br />mCircle = (MinimumBoundingCircle(geomColl)).getCircle()<br /><br /># Plot the geometries, BANG<br />plot(pntGeoms + [geomConvexH] + [mCircle])<br /></pre><br /><a onblur="try {parent.deselectBloggerImageGracefully();} catch(e) {}" href="http://1.bp.blogspot.com/_Z3xEawQVAXI/TCrecXTjfkI/AAAAAAAAAl0/InV-NJm7iU4/s1600/convex_mbc.png"><img style="display: block; margin: 0px auto 10px; text-align: center; cursor: pointer; width: 400px; height: 384px;" src="http://1.bp.blogspot.com/_Z3xEawQVAXI/TCrecXTjfkI/AAAAAAAAAl0/InV-NJm7iU4/s400/convex_mbc.png" alt="" id="BLOGGER_PHOTO_ID_5488443674778107458" border="0" /></a><br />I imported GeoScript modules and JTS Java classes in the example above. Jython can interoperate with Java libraries and classes because under the hood it’s Java. But how did I know what parameters a GeometryCollection accepted? I took a peek at the <a href="http://www.jarvana.com/jarvana/view/com/vividsolutions/jts/1.11/jts-1.11-javadoc.jar%21/com/vividsolutions/jts/geom/GeometryCollection.html">JTS Java doc</a> and followed the bouncing argument ball. The weaving of colorful geometries happens after that thanks to the great <a href="http://geoscriptblog.blogspot.com/2010/06/geometry-plotting.html">plotting functionality</a> from Justin.<br /><br />There are tricks and trapdoors to this unwrapped exploration - the syntax is not as clean and intuitive as working with GeoScript Python. But a lot of magic still happens. I try to use functions that return <a href="http://geoscript.org/learning/geom.html#learning-geom">geometry objects</a> so I can plot them with the geoscript.viewer.plot function. This way I can see each step visually. Another advantage is that JTS and GeoTools have very well organized Java docs (<a href="http://www.jarvana.com/jarvana/view/com/vividsolutions/jts/1.11/jts-1.11-javadoc.jar%21/overview-summary.html">here</a> and <a href="http://javadoc.geotools.fr/2.5/">there</a>). It’s hard to fail at finding your way with road maps this good.<br /><br />For example, mkgeomatics was deciding how to visually represent the output of a <a href="http://www.mkgeomatics.com/wordpress/">least-cost analysis</a> he performed on a Bellingham, WA dataset. Part of the cartographic workflow used Voronoi diagrams that were built in ArcGIS. Both of us wondered if we could leverage GeoScript Python and JTS to solve the Voronoi diagram work. After a little sandbox scripting I came up with the following test code. It outlines how to derive a Voronoi diagram from a Delaunay triangulation (A more succinct way to build Voronoi diagrams would<span style="font-weight: bold;"> </span>use the<span style="font-weight: bold;"> </span>com.vividsolutions.jts.triangulate.VoronoiDiagramBuilder class. But Dr JTS produced an <a href="http://lin-ear-th-inking.blogspot.com/2009/05/voronoi-diagrams-in-jts-111.html">informative post</a> on the relationship between both that I wanted to explore).<br /><pre class="brush:py">from geoscript.layer import Shapefile<br />from geoscript.viewer import plot<br />from com.vividsolutions.jts.geom import GeometryCollection<br />from com.vividsolutions.jts.geom import GeometryFactory<br />from com.vividsolutions.jts.geom import PrecisionModel<br />from com.vividsolutions.jts.algorithm import MinimumBoundingCircle<br />from com.vividsolutions.jts.triangulate import DelaunayTriangulationBuilder<br /><br />#Get a point shapefile<br />staLayer = Shapefile('/home/gregcorradini/wa_stations.shp')<br /><br />#Create a list of feature geometries<br />pntGeoms = [f.geom for f in staLayer.features()]<br /><br />#Create a GeometryCollection from our list of points<br />geomColl = GeometryCollection(pntGeoms,GeometryFactory(PreciscionModel(),4326))<br /><br />#Create a DelaunayTriangulationBuilder object and pass our geometry collection<br />dtb = DelaunayTriangulationBuilder()<br />dtb.setSites(geomColl)<br /><br />#Get the Delaunay triangle faces<br />triangles = dtb.getTriangles(GeometryFactory(PrecisionModel(),4326))<br /><br />#Make a pretty triangle picture<br />plot(pntGeoms + [triangles])<br /></pre><br /><a onblur="try {parent.deselectBloggerImageGracefully();} catch(e) {}" href="http://2.bp.blogspot.com/_Z3xEawQVAXI/TCrfi3suMUI/AAAAAAAAAl8/oalTy-Nnvd8/s1600/delaunayTriangles.png"><img style="display: block; margin: 0px auto 10px; text-align: center; cursor: pointer; width: 400px; height: 384px;" src="http://2.bp.blogspot.com/_Z3xEawQVAXI/TCrfi3suMUI/AAAAAAAAAl8/oalTy-Nnvd8/s400/delaunayTriangles.png" alt="" id="BLOGGER_PHOTO_ID_5488444886064443714" border="0" /></a><br /><pre class="brush:py">#Get the object's corresponding voronoi polygons in a java array<br />voronoiArray = dtb.getSubdivision().getVoronoiCellPolygons(GeometryFactory(PrecisionModel(),4326))<br /><br />#Make a pretty picture<br />plot(pntGeoms + [i for i in voronoiArray])<br /></pre><a onblur="try {parent.deselectBloggerImageGracefully();} catch(e) {}" href="http://1.bp.blogspot.com/_Z3xEawQVAXI/TCriLWm7XsI/AAAAAAAAAmE/iGAW5VogWGA/s1600/voronoi.png"><img style="display: block; margin: 0px auto 10px; text-align: center; cursor: pointer; width: 400px; height: 384px;" src="http://1.bp.blogspot.com/_Z3xEawQVAXI/TCriLWm7XsI/AAAAAAAAAmE/iGAW5VogWGA/s400/voronoi.png" alt="" id="BLOGGER_PHOTO_ID_5488447780579663554" border="0" /></a>That last Voronoi picture is hard to see in detail. So let's save off the Voronoi polygons to a shapefile and give each one the attributes of the original point that falls within it. The function 'changeGeometryType' is from a previous <a href="http://geoscriptblog.blogspot.com/2010/06/buffering-features-with-geoscript.html">post</a>.<br /><pre class="brush:py"><br />#Create a voronoi geom.Poly schema based on old layer schema<br />voronoiSchema = changeGeometryType('voronoiWAStations',geom.Polygon,staLayer)<br /><br />#Create voronoi layer<br />vLayer = ws.create(schema=voronoiSchema)<br /><br />#For each voronoi polygon test which point is within it and pass attributes<br />for poly in voronoiArray:<br /> for f in staLayer.features():<br /> <br /> #Is my point within the polygon?<br /> if f.geom.within(poly):<br /> <br /> #Copy all point feature attributes into a dictionary<br /> dictFeat = dict(f.iteritems())<br /> <br /> #Overwrite the_geom with voronoi polygon<br /> dictFeat['the_geom'] = poly<br /> <br /> #Create a new feature and add it to voronoi layer<br /> vLayer.add(voronoiSchema.feature(dictFeat))<br /></pre><br /><a onblur="try {parent.deselectBloggerImageGracefully();} catch(e) {}" href="http://2.bp.blogspot.com/_Z3xEawQVAXI/TC7B4R4YnoI/AAAAAAAAAmM/SsuFapdzfCw/s1600/final_voronoi.jpeg"><img style="display: block; margin: 0px auto 10px; text-align: center; cursor: pointer; width: 400px; height: 174px;" src="http://2.bp.blogspot.com/_Z3xEawQVAXI/TC7B4R4YnoI/AAAAAAAAAmM/SsuFapdzfCw/s400/final_voronoi.jpeg" alt="" id="BLOGGER_PHOTO_ID_5489538168427945602" border="0" /></a>The zoomed-in shot above shows the original points in pink, Delaunay triangles in gray and Voronoi diagram in green.Unknownnoreply@blogger.com0tag:blogger.com,1999:blog-8042770441502149016.post-39735292308022492512010-06-22T14:07:00.001-07:002010-06-25T12:38:06.244-07:00Geometry PlottingI often have found myself looking for a quick and easy way to visualize a geometry object. Maybe I have a handful of geometry objects and I want to visualize some relationship (for example intersection) that exists between them. <br /><br />Recently added to geoscript python is the ability to create an xy plot from a set of geometry objects. This functionality comes courtesy of the <a href="http://github.com/jdeolive/geoscript-py/blob/master/geoscript/viewer.py#L78">plot</a> function that is located in the viewer module.<br /><br /><pre class="brush:py"><br />from geoscript.viewer import plot<br />from geoscript import geom<br /><br />g = geom.LineString((0,0), (10,10))<br />plot(g)<br /></pre><br /><br /><img src="http://img.skitch.com/20100622-dxq3hg9ec83c9gcrdd8kymihwq.jpg" /><br /><br />Nothing all that earth shattering but something can be utilized in a variety of different ways. For instance visualizing an intersection as mentioned above:<br /><br /><pre class="brush:py"><br />g1 = geom.readWKT('POLYGON((-87.1875 49.5703125, -124.453125 -29.8828125, -92.8125 -58.7109375, -43.59375 -69.9609375, -4.21875 -60.1171875, 16.171875 -29.8828125, 16.171875 -5.2734375, 0 27.7734375, -20.390625 48.8671875, -55.546875 53.7890625, -78.046875 54.4921875, -87.1875 49.5703125))')<br /><br />g2 = geom.readWKT('POLYGON((24.609375 51.6796875, -28.125 24.9609375, -25.3125 -20.7421875, 24.609375 -61.5234375, 74.53125 -51.6796875, 108.984375 -23.5546875, 92.109375 27.0703125, 61.171875 55.1953125, 40.078125 55.1953125, 24.609375 51.6796875))')<br /><br />plot([g1.intersection(g2), g1,g2])<br /></pre><br /><br /><img src="http://img.skitch.com/20100622-en416rwwegpuutnig1ine2s2ka.jpg"/><br /><br />The plotting functionality makes use of the <a href="http://www.jfree.org/jfreechart/"/>JFreeChart</a> library, a popular open source Java framework for creating diagrams and charts. Another example of one of the benefits of the marriage of Java and Python that is Jython.<br /><br />A recent use I made of this new functionality was with regard to geometry simplification. I wanted to quickly visualize a simplified geometry at different distance tolerances. Using everybody's favourite layer as an example:<br /><br /><pre class="brush:py"><br />from geoscript.layer import Shapefile<br />from geoscript.geom import simplify<br /><br />shp = Shapefile('tests/work/states.shp')<br /><br />texas = [f for f in shp.features("STATE_NAME = 'Texas'")][0].geom<br />plot(texas)<br /><br />texas_simple = simplify(texas, 0.1)<br />plot(texas_simple)<br /><br />texas_simple2 = simplify(texas, 0.5)<br />plot(texas_simple2)<br /></pre><br /><br /><img src="http://img.skitch.com/20100622-bmm34r9raun3f34s5r8xee6pq1.jpg"/><br /><br />The simplification routine is the Douglas-Peucker algorithm provided out of the box by JTS.Justin Deoliveirahttp://www.blogger.com/profile/04340292357363191337noreply@blogger.com2tag:blogger.com,1999:blog-8042770441502149016.post-72367148772910831782010-06-20T15:14:00.000-07:002010-06-27T06:31:08.123-07:00Saving Weather Stations to a Shapefile with PythonThe National Weather Service offers several XML <a href="http://www.nws.noaa.gov/xml/current_obs/">data feeds</a> related to weather conditions. In this tutorial we will use the <a href="http://www.weather.gov/xml/current_obs/index.xml">list</a> of all national observation stations to create a shapefile of Washington stations with Python.<br /><pre class="brush:py"><br />def createXMLfileFromURL(strURL,strOutFilePath):<br /><br /> #Open url and read<br /> res = urllib2.urlopen(strURL)<br /> data = res.read()<br /> <br /> #Write data to file<br /> writer = open(strOutFilePath,'w')<br /> writer.write(data)<br /> writer.close()<br /><br />if __name__=='__main__':<br /><br /> #Grab the XML feed and write it to a file<br /> noaaIndexXML = '/home/gregcorradini/XMLfeeds/indexNOAAII.xml'<br /> createXMLfileFromURL('http://www.weather.gov/xml/current_obs/index.XML',noaaIndexXML)<br /><br /> #Use ElementTree to walk XML for 'station'<br /> root= (etree.parse(noaaIndexXML)).getroot()<br /> listStations = root.findall('station')<br /><br /> #Findall WA stations in XML document<br /> waStations = [i for i in listStations if i.findall('state')[0].text == 'WA']<br /><br /> #Create a dictionary for each station<br /> dictWa = {}<br /> for i in waStations:<br /> #Notice we are using the lat/lng values to create a geom.Point<br /> dictWa[i.findall('station_id')[0].text] = {'the_geom':geom.Point(float(i.findall('longitude')[0].text),float(i.findall('latitude')[0].text)),'station_id':i.findall('station_id')[0].text,'lat':float(i.findall('latitude')[0].text), 'lng':float(i.findall('longitude')[0].text), 'xml_url': i.findall('xml_url')[0].text}<br /><br /> #Create a schema for weather stations<br /> staSchema = schema.Schema('wa_stations_notemp',[('the_geom',geom.Point,Projection('epsg:4326')),('station_id',str),('lat',float),('lng',float),('xml_url',str)])<br /><br /> #Get a workspace<br /> ws = Directory('/home/gregcorradini/GeoTools/geoscript_wrk/isolines/data/US_shapefiles/')<br /><br /> #Create a stations shapefile with our schema<br /> staLayer = ws.create('wa_stations_notemp',schema=staSchema)<br /><br /> #Create a feature for each station value<br /> counter = 1<br /> for key in dictWa.keys():<br /> staFeature = staLayer.schema.feature(dictWa[key],str(counter))<br /> staLayer.add(staFeature)<br /> counter += counter<br /></pre><br /><a onblur="try {parent.deselectBloggerImageGracefully();} catch(e) {}" href="http://4.bp.blogspot.com/_Z3xEawQVAXI/TB6Yfm0hU9I/AAAAAAAAAlk/ykl2BAZGw3c/s1600/xml2shp.jpeg"><img style="display: block; margin: 0px auto 10px; text-align: center; cursor: pointer; width: 374px; height: 263px;" src="http://4.bp.blogspot.com/_Z3xEawQVAXI/TB6Yfm0hU9I/AAAAAAAAAlk/ykl2BAZGw3c/s400/xml2shp.jpeg" alt="" id="BLOGGER_PHOTO_ID_5484989064948044754" border="0" /></a>Unknownnoreply@blogger.com2tag:blogger.com,1999:blog-8042770441502149016.post-68635161591141447222010-06-15T04:40:00.000-07:002010-06-15T05:31:26.243-07:00Merging Shapefiles with JavaScript<a href="http://geoscript.org">GeoScript</a> aims to provide a framework for scripting common data processing tasks with <a href="http://geotools.org">GeoTools</a> and the <a href="http://tsusiatsoftware.net/jts/main.html">Java Topology Suite</a>. Much in the same way that scripting languages like Python can leverage the powerful OGR/GDAL libraries, GeoScript offers similar functionality from languages that run on the JVM (currently <a href="http://geoscript.org/groovy/index.html">Groovy</a>, <a href="http://geoscript.org/js/index.html">JavaScript</a>, <a href="http://geoscript.org/py/index.html">Python</a>, or <a href="http://geoscript.org/scala/index.html">Scala</a>).<br /><br />In a recent post on his blog, <a href="http://darrencope.com/2010/05/07/merge-a-directory-of-shapefiles-using-ogr/">Darren Cope wrote</a> about merging a directory of shapefiles, a common data processing task. Here I'll demonstrate how to do the same thing with the JavaScript version of GeoScript.<br /><br />This example uses the <a href="http://geoscript.org/js/api/workspace.html">workspace</a> and <a href="http://geoscript.org/js/api/layer.html">layer</a> modules. I'm working with a directory containing shapefiles of US Census block groups for each state, merging them into a single shapefile for the country.<br /><br /><pre class="brush:javascript; gutter:false"><br />// import modules<br />var workspace = require("geoscript/workspace");<br />var layer = require("geoscript/layer");<br /><br />// create workspaces from existing directories<br />var source = new workspace.Directory("path/to/source_dir");<br />var target = new workspace.Directory("path/to/target_dir");<br /><br />// iterate through layers in source workspace<br />var country;<br />source.names.forEach(function(name) {<br /> // create state layer from existing shapefile<br /> var state = source.get(name);<br /> // create country layer first time through<br /> if (!country) {<br /> country = new layer.Layer({<br /> schema: state.schema.clone({name: "country"})<br /> });<br /> // this creates the new shapefile on disk<br /> target.add(country);<br /> }<br /> // iterate through source features to add each to target <br /> state.features.forEach(function(feature) {<br /> country.add(feature);<br /> });<br />});<br /></pre><br /><br />While it is hard to beat the terseness of the <a href="http://darrencope.com/2010/05/07/merge-a-directory-of-shapefiles-using-ogr/">ogr2ogr example</a>, the JavaScript version of GeoScript provides syntax that I hope is familiar to a <a href="http://www.cio.com/article/446829/PHP_JavaScript_Ruby_Perl_Python_and_Tcl_Today_The_State_of_the_Scripting_Universe?contentId=446829&slug=&">larger population</a> of scripters.Tim Schaubhttp://www.blogger.com/profile/12243618766226981189noreply@blogger.com0tag:blogger.com,1999:blog-8042770441502149016.post-17617279086513267372010-06-14T08:39:00.000-07:002010-06-14T20:39:45.045-07:00Calculating Convex Hull and Minimum Bounding Circle with GroovyThe Java Topology Suite, which <a href="http://geoscript.org/">Groovy GeoScript</a> wraps, contains spatial operators that act on a group of Features or Geometries. In this post, I collect all geometries from a shapefile to calculate the convex hull and minimum bounding circle. This post builds on previous blog entries where I use GeoScript to extract <a href="http://geoscriptblog.blogspot.com/2010/04/calculating-centroids-with-geoscript.html">centroids</a> and <a href="http://geoscriptblog.blogspot.com/2010/05/buffering-features-with-groovy.html">buffer</a> Features.<br /><br /><b>Convex Hull</b><br /><pre class="brush:groovy">// Import GeoScript modules<br />import geoscript.layer.*<br />import geoscript.feature.*<br />import geoscript.geom.*<br /><br />// Get the shapefile<br />Shapefile shp = new Shapefile('states_centroids.shp');<br /><br />// Create a new Schema<br />Schema schema = new Schema('states_convex_hull', [['the_geom','Polygon','EPSG:4326']])<br /><br />// Create our new Layer<br />Layer layer = shp.workspace.create(schema)<br /><br />// Collect the Geometries<br />List geoms = shp.features.collect{f->f.geom}<br /><br />// Create a GeometryCollection from the List of Geometries<br />GeometryCollection geomCol = new GeometryCollection(geoms)<br /><br />// Get the Convex Hull from the GeometryCollection<br />Geometry convexHullGeom = geomCol.convexHull<br /><br />// Add the Convex Hull Geometry as a Feature<br />layer.add(schema.feature([convexHullGeom]))<br /></pre><br /><a onblur="try {parent.deselectBloggerImageGracefully();} catch(e) {}" href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjZtiIq4E36DhxE3gZPbSKsrBQG1vbJL5KTo7lKlro71ShcH6ns2HIM8ZTijLcDQHve5TwJ5SlcBOd2rnIUca0fVVhOVwKlz9Qyo3wcCpbEXciAnqZSTW-jsAGx_8wRY1XY1_kZsn8OhZ81/s1600/convex_hull.png"><img style="display:block; margin:0px auto 10px; text-align:center;cursor:pointer; cursor:hand;width: 400px; height: 180px;" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjZtiIq4E36DhxE3gZPbSKsrBQG1vbJL5KTo7lKlro71ShcH6ns2HIM8ZTijLcDQHve5TwJ5SlcBOd2rnIUca0fVVhOVwKlz9Qyo3wcCpbEXciAnqZSTW-jsAGx_8wRY1XY1_kZsn8OhZ81/s400/convex_hull.png" border="0" alt="" id="BLOGGER_PHOTO_ID_5475985887301440946"></a><br /><br /><b>Minimum Bounding Circle</b><br /><pre class="brush:groovy"><br />// Import GeoScript modules<br />import geoscript.layer.*<br />import geoscript.feature.*<br />import geoscript.geom.*<br /><br />// Get the Shapefile<br />Shapefile shp = new Shapefile('states_centroids.shp')<br /><br />// Create a new Schema<br />Schema schema = new Schema('states_minimum_bounding_circle', [['the_geom','Polygon','EPSG:4326']])<br /><br />// Create the new Layer<br />Layer layer = shp.workspace.create(schema)<br /><br />// Collect Geometries from the Shapefile<br />List geoms = shp.features.collect{f->f.geom}<br /><br />// Create a GeometryCollection from the List of Geometries<br />GeometryCollection geomCol = new GeometryCollection(geoms)<br /><br />// Get the Minimum Bounding Circle from the GeometryCollection<br />Geometry circleGeom = geomCol.minimumBoundingCircle<br /><br />// Add the Minimum Bounding Circle Geometry as a Feature<br />layer.add(schema.feature([circleGeom]))<br /></pre><br /><a onblur="try {parent.deselectBloggerImageGracefully();} catch(e) {}" href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiXIVrk2IUM_7wthnZLCtc_nYPxGkotKV9PF6KrXkOgtkxxLrbECvrTVncaUJQKet7X0vIz7yUWv9OGTJdjhpk1GsmNdgfdE8SE86wy4XHBQvPP9Lg0bS2QJirDRsXy30il6A4nwEYLmoRD/s1600/states_minboundingcircle.png"><img style="display:block; margin:0px auto 10px; text-align:center;cursor:pointer; cursor:hand;width: 400px; height: 331px;" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiXIVrk2IUM_7wthnZLCtc_nYPxGkotKV9PF6KrXkOgtkxxLrbECvrTVncaUJQKet7X0vIz7yUWv9OGTJdjhpk1GsmNdgfdE8SE86wy4XHBQvPP9Lg0bS2QJirDRsXy30il6A4nwEYLmoRD/s400/states_minboundingcircle.png" border="0" alt="" id="BLOGGER_PHOTO_ID_5475986007467915522"></a>Jared Ericksonhttp://www.blogger.com/profile/13350047301477401257noreply@blogger.com2tag:blogger.com,1999:blog-8042770441502149016.post-7536707011624578402010-06-12T13:59:00.000-07:002010-06-26T09:08:42.949-07:00Buffering Features with GeoScript PythonIn a previous <a href="http://geoscriptblog.blogspot.com/2010/05/buffering-features-with-groovy.html">post</a>, Jared buffered points and added them to a polygon shapefile. Here's the workflow in Python GeoScript.<br /><pre class="brush:py"><br />from geoscript import proj<br />from geoscript import geom<br />from geoscript.layer import Shapefile<br />from geoscript.feature import schema<br /><br />def changeGeometryType(strName,geomType,shpParent):<br /><br /> #Get all shapefile schema attributes except 'the_geom'<br /> parentS = [(str(i.name),i.typ) for i in shpParent.schema.fields if not issubclass(i.typ,geom.Geometry)]<br /><br /> #Get shapefile Projection<br /> pPrj = shpParent.schema.get('the_geom').proj<br /><br /> #Reinsert geomType into schema<br /> parentS.insert(0,('the_geom',geomType,pPrj))<br /> <br /> #Create centroid buffer schema<br /> centroidSchema = schema.Schema(strName,parentS)<br /><br /> return centroidSchema<br /><br /><br />if __name__=='__main__':<br /><br /> #Get centroid shapefile<br /> shp = Shapefile('/home/gregcorradini/GeoTools/geoscript_wrk/centroid_buffer/data/centroids.shp')<br /><br /> #Create a geom.Polygon schema from centroid schema<br /> centroidSchema = changeGeometryType('centroid_buffer',geom.Polygon,shp)<br /><br /> #Get shapefile workspace<br /> ws = shp.workspace<br /><br /> #Create a buffered centroid layer based on centroid schema<br /> buffCentLayer = ws.create(schema=centroidSchema)<br /><br /> #For each centroid feature, create a buffered centroid, add it to buffered layer<br /> for f in shp.features():<br /> <br /> #Copy all feature attributes into a dictionary<br /> dictFeat = dict(f.iteritems())<br /> <br /> #Buffer 'the_geom' 3 degrees<br /> dictFeat['the_geom'] = f.geom.buffer(3)<br /><br /> print dictFeat<br /> <br /> #Create a buffered centroid feature with new dictionary attributes<br /> buffFeature = buffCentLayer.schema.feature(dictFeat,f.id)<br /><br /> #Add new feature to buffered layer<br /> buffCentLayer.add(buffFeature)<br /></pre><br /><br /><a onblur="try {parent.deselectBloggerImageGracefully();} catch(e) {}" href="http://2.bp.blogspot.com/_Z3xEawQVAXI/TB4rMJqtdRI/AAAAAAAAAlU/X7TbcPXpus4/s1600/centroidBuffers.jpeg"><img style="margin: 0px auto 10px; display: block; text-align: center; cursor: pointer; width: 374px; height: 299px;" src="http://2.bp.blogspot.com/_Z3xEawQVAXI/TB4rMJqtdRI/AAAAAAAAAlU/X7TbcPXpus4/s400/centroidBuffers.jpeg" alt="" id="BLOGGER_PHOTO_ID_5484868883937326354" border="0" /></a>Unknownnoreply@blogger.com0tag:blogger.com,1999:blog-8042770441502149016.post-80891334200674887582010-06-12T11:18:00.000-07:002010-06-26T08:54:45.343-07:00Calculating Centroids with GeoScript Python<a href="http://geoscriptblog.blogspot.com/2010/04/calculating-centroids-with-geoscript.html">Jared</a> and <a href="http://geoscriptblog.blogspot.com/2010/05/calculating-centroids-with-javascript.html">Tim</a> have shown us how to create centroids from a polygon shapefile in two different GeoScript implementations. Now let's see what the code looks like in Python GeoScript.<br /><pre class="brush:py"><br />from geoscript import proj<br />from geoscript import geom<br />from geoscript.layer import Shapefile<br />from geoscript.feature import schema<br /><br />def changeGeometryType(strName,geomType,shpParent):<br /><br /> #Put shapefile schema attributes except 'the geom' in a list of tuples<br /> parentS = [(str(i.name),i.typ) for i in shpParent.schema.fields if not issubclass(i.typ,geom.Geometry)]<br /><br /> #Get shapefile Projection<br /> pPrj = shpParent.schema.field('the_geom').proj<br /><br /> #reinsert geomType and proj into schema<br /> parentS.insert(0,('the_geom',geomType,pPrj))<br /> <br /> #create centroid schema with the list<br /> centroidSchema = schema.Schema(strName,parentS)<br /><br /> return centroidSchema<br /><br /><br />if __name__=='__main__':<br /><br /> #Get states shapefile<br /> shp = Shapefile('/home/gregcorradini/GeoTools/geoscript_wrk/centroid_buffer/data/US_shapefiles.shp')<br /><br /> #Substitute geom.Point into schema<br /> centroidSchema = changeGeometryType('centroids',geom.Point,shp)<br /><br /> #Get workspace<br /> ws = shp.workspace<br /><br /> #Create a centroid layer based on centroidSchema<br /> centroidLayer = ws.create(schema=centroidSchema)<br /><br /> #Get all state features<br /> stateFeatures = [f for f in shp.features()]<br /><br /> #For each state polygon get its centroid and add it to our centroid layer<br /> for f in stateFeatures:<br /> <br /> #create dictionary to hold each feature's attributes<br /> newAtts = {}<br /> <br /> for key in f.attributes.keys():<br /> if issubclass(f.attributes[key].__class__,geom.Geometry):<br /> newAtts[str(key)] = f.geom.centroid<br /> else:<br /> newAtts[str(key)] = f.attributes[key]<br /><br /> print newAtts<br /><br /> #Create a feature based on dictionary<br /> centFeature = centroidLayer.schema.feature(newAtts,f.id)<br /><br /> #Add it to centroid layer<br /> centroidLayer.add(centFeature)<br /><br /><br /></pre><br /><br /><a onblur="try {parent.deselectBloggerImageGracefully();} catch(e) {}" href="http://1.bp.blogspot.com/_Z3xEawQVAXI/TB4n7d0H22I/AAAAAAAAAlE/kmSy8qaWJAc/s1600/centroids.jpeg"><img style="margin: 0px auto 10px; display: block; text-align: center; cursor: pointer; width: 374px; height: 299px;" src="http://1.bp.blogspot.com/_Z3xEawQVAXI/TB4n7d0H22I/AAAAAAAAAlE/kmSy8qaWJAc/s400/centroids.jpeg" alt="" id="BLOGGER_PHOTO_ID_5484865298752854882" border="0" /></a>Unknownnoreply@blogger.com0tag:blogger.com,1999:blog-8042770441502149016.post-67303726363331880472010-05-10T11:27:00.000-07:002010-05-10T12:03:03.494-07:00Calculating Centroids with JavaScript<a onblur="try {parent.deselectBloggerImageGracefully();} catch(e) {}" href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEj478tcr98lMHKlLHriSOPANGL1hX-llw3td7b9Ay6-K5k8OIGlXvDzmORPYRPBoNXJA5Apt8TaNBRxQhuzYhSwpsO6ujGrf1X9ZozaFScUjoF6dbEW6QUQEVl8KBP8NZnpGx1CVpzHbwg/s1600/centroids.gif"><img style="display:block; margin:0px auto 10px; text-align:center;cursor:pointer; cursor:hand;width: 320px; height: 159px;" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEj478tcr98lMHKlLHriSOPANGL1hX-llw3td7b9Ay6-K5k8OIGlXvDzmORPYRPBoNXJA5Apt8TaNBRxQhuzYhSwpsO6ujGrf1X9ZozaFScUjoF6dbEW6QUQEVl8KBP8NZnpGx1CVpzHbwg/s320/centroids.gif" border="0" alt=""id="BLOGGER_PHOTO_ID_5469716759809437074" /></a><br /><br />In a <a href="http://geoscriptblog.blogspot.com/2010/04/calculating-centroids-with-geoscript.html">previous post</a>, Jared described how to use the <a href="http://geoscript.org/groovy/index.html">Groovy GeoScript</a> to create a new Shapefile based on centroids of features in an existing Shapefile. Today, I'll do the same using the <a href="http://geoscript.org/js/index.html">JavaScript</a> implementation.<br /><br />One significant difference from the other implementations is that there is only one layer type in the JavaScript flavor of GeoScript. To create a layer representing a Shapefile, I'll create a directory workspace and get the layer from that.<br /><br />In this example, I'll also use the <code>clone</code> methods to create a new layer schema and new features. These <code>clone</code> methods allow for new objects to be created with modified properties. This is useful in creating the schema for the centroids layer, because we want the same fields as the states layer with the exception of the geometry field. In this case, we're creating a layer with point geometries instead of multi-polygon geometries.<br /><br />Without further ado, here is the code:<br /><pre class="brush:javascript; gutter:false"><br />// import Directory and Layer constructors<br />var Directory = require("geoscript/workspace").Directory;<br />var Layer = require("geoscript/layer").Layer;<br /><br />// create a directory workspace from an existing directory on disk<br />var dir = new Directory("data/shapefiles");<br /><br />// create a layer based on an existing shapefile in the directory<br />var states = dir.get("states");<br /><br />// create a new schema with a Point geometry type instead of MultiPolygon<br />var schema = states.schema.clone({<br /> // give the schema a new name<br /> name: "centroids", <br /> fields: [<br /> // overwrite existing field named "the_geom"<br /> {name: "the_geom", type: "Point", projection: "EPSG:4326"}<br /> ]<br />});<br /><br />// create a new temporary layer with the new schema<br />var centroids = new Layer({<br /> schema: schema<br />});<br /><br />// add the layer to existing workspace (this creates a new shapefile on disk)<br />dir.add(centroids);<br /><br />// iterate through the state features to create features with state centroids<br />states.features.forEach(function(state) {<br /> var centroid = state.clone({<br /> schema: schema,<br /> values: {<br /> // overwrite the geometry value with the state centroid<br /> the_geom: state.geometry.centroid<br /> }<br /> });<br /> // add the new feature to the layer (this writes to the shapefile)<br /> centroids.add(centroid);<br />});<br /></pre><br /><br />The GeoScript <code>viewer</code> module exports a <code>draw</code> method. You can use this to render a layer, or arrays of geometies or features.<br /><br /><pre class="brush:javascript; gutter:false"><br />// use the viewer module to draw collections of features or geometries<br />var draw = require("geoscript/viewer").draw;<br />var geometries = [];<br /><br />// add all state geometries to the array<br />states.features.forEach(function(state) {<br /> geometries.push(state.geometry);<br />});<br /><br />// add buffered centroids so they will be visible<br />centroids.features.forEach(function(centroid) {<br /> geometries.push(centroid.geometry.buffer(0.5));<br />});<br /><br />// draw all geometries<br />draw(geometries);<br /></pre><br /><br />You'll find this example in the <a href="http://github.com/tschaub/geoscript-js/tree/master/examples/">examples directory</a> of the repository. See the <a href="http://geoscript.org/js/api/index.html">JavaScript API documentation</a> for more detail on using GeoScript.Tim Schaubhttp://www.blogger.com/profile/12243618766226981189noreply@blogger.com0tag:blogger.com,1999:blog-8042770441502149016.post-28672080998844954682010-05-07T19:44:00.000-07:002010-05-10T21:28:56.678-07:00Buffering Features with GroovyIn my <a href="http://geoscriptblog.blogspot.com/2010/04/calculating-centroids-with-geoscript.html">last post</a>, I calculated centroids from one shapefile and saved them to another using a <a href="http://geoscript.org/">GeoScript</a> <a href="http://groovy.codehaus.org/">Groovy</a> script. This time I buffer these centroids and save them to a polygon shapefile. Since GeoScript is based on the <a href="http://sourceforge.net/projects/jts-topo-suite/">Java Topology Suite</a> library you can take advantage of any of its geometry operations - intersection, union and difference.<br /><pre class="brush:groovy">// Import Geoscript modules<br />import geoscript.layer.*<br />import geoscript.feature.*<br />import geoscript.geom.*<br /><br />// Get the Shapefile we want to buffer<br />Shapefile shp = new Shapefile('states_centroids.shp')<br /><br />// Create a new Schema but with Polygon Geometry<br />Schema schema = shp.schema.changeGeometryType('Polygon','states_buffers')<br /><br />// Create the new Layer<br />Layer bufferLayer = shp.workspace.create(schema)<br /><br />// Specify the buffer distance <br />double distance = 2 // decimal degrees<br /><br />// Iterate through each Feature using a closure<br />shp.features.each{f -><br /><br /> // Create a Map for the new attributes<br /> Map attributes = [:]<br /><br /> // Set attribute values<br /> f.attributes.each{k,v -><br /> if (v instanceof Geometry) {<br /> attributes[k] = v.buffer(distance)<br /> }<br /> else {<br /> attributes[k] = v<br /> }<br /> }<br /><br /> // Create a new Feature with the new attributes<br /> Feature feature = schema.feature(attributes, f.id)<br /><br /> // Add it to the buffer Layer<br /> bufferLayer.add(feature)<br />}<br /></pre><br />The Layer.getCursor() method, used in the last post, reads one Feature at a time. In this post I use Layer.getFeatures(). This method reads all of the Features in to a list. Here is what you get:<br /><br /><a onblur="try {parent.deselectBloggerImageGracefully();} catch(e) {}" href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgmFXI4aAQ5H_6VxJBA_Acr_BEKTpFEGvDynFjlXyDL2j4QR8hpKerVIP81_nXtc30K-gr52P12x6i_iAqo7XOMTlAm5D9H5u9YRXcwp7zn-z1dCF74vC4grZvbOIbNuM1nRSEOvyRJ02rP/s1600/states_buffer.png"><img style="display:block; margin:0px auto 10px; text-align:center;cursor:pointer; cursor:hand;width: 400px; height: 179px;" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgmFXI4aAQ5H_6VxJBA_Acr_BEKTpFEGvDynFjlXyDL2j4QR8hpKerVIP81_nXtc30K-gr52P12x6i_iAqo7XOMTlAm5D9H5u9YRXcwp7zn-z1dCF74vC4grZvbOIbNuM1nRSEOvyRJ02rP/s400/states_buffer.png" border="0" alt="" id="BLOGGER_PHOTO_ID_5468725151105414210"></a>Jared Ericksonhttp://www.blogger.com/profile/13350047301477401257noreply@blogger.com0tag:blogger.com,1999:blog-8042770441502149016.post-71503550729312559932010-04-29T19:52:00.000-07:002010-05-04T20:47:11.062-07:00Calculating Centroids with GeoScript Groovy<div style="text-align: left;"><span class="Apple-style-span" style=" ;font-size:medium;"><a href="http://geoscript.org/">GeoScript</a> uses the <a href="http://www.geotools.org/">GeoTools</a> and <a href="http://sourceforge.net/projects/jts-topo-suite/">Java Topology Suite</a> libraries to provide easy-to-use scripting APIs with implementations in <a href="http://github.com/jdeolive/geoscript-py">P</a><a href="http://github.com/jdeolive/geoscript-py">ython</a>, <a href="http://github.com/tschaub/geoscript-js">JavaScript</a>, <a href="http://github.com/dwins/geoscript.scala">Scala</a> and <a href="http://github.com/jericks/geoscript-groovy">Groovy</a>. The GeoScript web site provides simple code snippets in all four languages to help you get started. In this post, I address a real world example: calculating centroids from one GIS layer and saving them to another layer using Groovy. This example highlights the GeoScript Layer, Feature, and Geometry modules.</span></div><pre class="brush:groovy">// import GeoScript modules<br />import geoscript.layer.*<br />import geoscript.feature.*<br />import geoscript.geom.*<br /><br />// This is the Shapefile we want to extract centroids from<br />Shapefile shp = new Shapefile('states.shp')<br /><br />// Create a new Schema (based on the above Shapefile, but with Point Geometry)<br />Schema schema = shp.schema.changeGeometryType('Point','states_centroids')<br /><br />// Create the new centroid Layer<br />Layer centroidLayer = shp.workspace.create(schema)<br /><br />// Use the Cursor to loop through each Feature in the Shapefile<br />Cursor cursor = shp.cursor<br />while(cursor.hasNext()) {<br /><br /> // Get the next Feature<br /> Feature f = cursor.next()<br /> <br /> // Create a Map for the new attributes<br /> Map attributes = [:]<br /><br /> // For each attribute in the shapefile<br /> f.attributes.each{k,v -><br /><br /> // If its Geometry, find the centroid<br /> if (v instanceof Geometry) {<br /> attributes[k] = v.centroid<br /> }<br /> // Else, they stay the same<br /> else {<br /> attributes[k] = v<br /> }<br /> }<br /><br /> // Create a new Feature with the new attributes<br /> Feature feature = schema.feature(attributes, f.id)<br /><br /> // Add it to the centroid Layer<br /> centroidLayer.add(feature)<br />}<br /><br />// Always remember to close the cursor<br />cursor.close()<br /></pre><div>Here is what you get when you run the script against a shapefile of the 48 contiguous states (sorry Alaska and Hawaii):</div><br /><a onblur="try {parent.deselectBloggerImageGracefully();} catch(e) {}" href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhaOFgfgOwmpBqfK4fQrzWDErAqnqTFrDBoUS5X5QB0ZyKmvgm0b6bGovHtFuk8gv-SooebZFreAgRacvoy51xwMlK638pzxSSfyig9V_WhgJ_I4nxQ2ir3p_Q_jtEunnqITu42788YtV8t/s1600/centroids_states.png" style="text-decoration: none;"><img style="display:block; margin:0px auto 10px; text-align:center;cursor:pointer; cursor:hand;width: 400px; height: 164px;" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhaOFgfgOwmpBqfK4fQrzWDErAqnqTFrDBoUS5X5QB0ZyKmvgm0b6bGovHtFuk8gv-SooebZFreAgRacvoy51xwMlK638pzxSSfyig9V_WhgJ_I4nxQ2ir3p_Q_jtEunnqITu42788YtV8t/s400/centroids_states.png" border="0" alt="" id="BLOGGER_PHOTO_ID_5465764888433611890" /></a>Jared Ericksonhttp://www.blogger.com/profile/13350047301477401257noreply@blogger.com0tag:blogger.com,1999:blog-8042770441502149016.post-14060285536437001642010-03-26T07:52:00.000-07:002010-03-26T09:16:03.313-07:00WelcomeThe GeoScript team is happy to launch the official blog for the project. This will be the place to watch for all news and announcements related to releases, new features, and general project happenings. Welcome!Justin Deoliveirahttp://www.blogger.com/profile/04340292357363191337noreply@blogger.com0