TL;DR: Generating stable RPCs is hard, and I'm not sure how to do it!
RPCs in the context of the geospatial world are equations for describing how a geometrically raw sensor image can be related to a location in the real world. RPC can variously expand to Rational Polynomial Coefficient or Rapid Positioning Coefficients. They consist of 80 coefficents for a ratio of two third degree polynomials (for x and y) plus 10 more values for normalizing image and geographic coordinates.
A couple background papers include one by Vincent Tao and others on Understanding the Rational Function Model (pdf), and one by Gene Dial and J. Grodecki titled RPC Replacement Camera Models (pdf). These date from the mid-2000s, but they reference some earlier work by Gene Dial and J. Grodecki who seem to have played the seminal role in promoting the use of RPCs as generic replacements for rigerous camera models. I believe the major satellite image vendors have now been delivering RPC models (example) for raw scenes for on the order of a decade.
I, and others, have added support to GDAL for reading and representing RPCs as metadata over the a number of years, as well as support for evaluating RPCs in gdalwarp as a way of georeferencing images. At it's very simplest a command like the following can be used to roughly rectify a raw image that has RPCs.
gdalwarp -rpc raw.tif rectified.tif
In fact, at FOSS4G NA this year, a command like this was the big reveal in Paul Morin's discussion of how the Polar Geospatial Center was able to efficiently process huge numbers of satellite scenes of the polar regions. I was rather surprised since I wasn't aware of this capability being in significant use. However the above use example will usually do a relatively poor job. RPCs encapsulate the perspective view of the satellite and as such you need to know the elevation of an early location in order to figure out exactly where it falls on the raw satellite image. If you had a relatively flat image at a known elevation you can specify a fixed elevation like this:
gdalwarp -rpc -to RPC_HEIGHT=257 raw.tif rectified.tif
However, in an area that is not so flat, you really need a reasonably could DEM (digital elevation model) of the region so that distortion effects from terrain can be taken into account. This is particularly important if the view point is relatively low (ie. for an airphoto) or if the image is oblique rather than nadir (straight down). You should be able to do this using a command referencing a DEM file for the area like this, though to be honest I haven't tried it recently (ever?) and I'm not sure how well it works:
gdalwarp -rpc -to RPC_DEM=srtm.tif raw.tif rectified.tif
The above is essentially review from my point of view. But with my recent move from Google to Planet Labs I have also become interested in how to generate RPCs given knowledge of a camera model, and the perspective view from a satellite. This would give us the opportunity to deliver raw satellite scenes with a generic description of the geometry that can already be exploited in a variety of software packages, including GDAL, OSSIM, Orfeo Toolbox and proprietary packages like ArcGIS, Imagine and PCI Geomatica. However, I had no idea how this would be accomplished and while I earned an Honours BMath, I have let the modest amount of math I mastered in the late 80's rot ever since.
So, using what I am good at, I started digging around the web for open source implementations for generating RPCs. Unfortunately, RPC is a very generic term much more frequently used to mean "Remote Procedure Call" in the software development world. Nevertheless, after digging I learned that the Orfeo Toolbox (OTB) has a program for generating RPCs from a set of GCPs (ground control points). Ground control points are just a list of locations on an image where the real world location is known. That is a list of information like (pixel x, pixel y, longitude, latitude, height). That was exciting, and a bit of digging uncovered the fact that in fact Orfeo Toolbox's RPC capability is really just a wrapper for services used from OSSIM. In particular core RPC "solver" is found in the source files ossimRpcSolver.cc, and ossimRpcSolver.h. They appear to have been authored by the very productive Garrett Potts.
Based on the example of the Orfeo Toolbox code, I wrote a gcps2rpc C++ program that called out to libossim to transform GCPs read from a text file into RPCs written to a GDAL VRT file.
For my first attempt I created 49 GCPs in a 7x7 pattern. The elevations were all zero (ie. they were on the WGS84 ellipsoids) and they were generated using custom code mapping from satellite image pixels to ground locations on the ellipsoid. I included simple code in gcps2rpc to do a forward and reverse evaluation of the ground control points to see how well they fit the RPC. In the case of this initial experiment the forward transformation, (long, lat, height) to (pixel, line) worked well with errors in the order of 1/1000th of a pixel. Errors in the inverse transformation, (pixel, line, height) to (long, lat) was less good with worst case errors on the order of 1/10th of a pixel.
At first I wasn't sure what to make of the errors in the inverse direction. 1/10th of a pixel wasn't terrible from my perspective, but it was still a fairly dramatic difference from the tiny errors in the forward direction. I tried comparing the results from evaluating the RPCs with OSSIM and with GDAL. To test with GDAL I just ran gdaltransform something like this.
gdaltransform -rpc -i simulated_rgb.vrt
-85.0495014028 39.7243133317 (entered in terminal)
2016.00011619623 1792.00003975369 0 (output in terminal)
and:
gdaltransform -rpc simulated_rgb.vrt
2016 1792 (entered in terminal)
-85.0495014329239 39.7243132735733 0 (output in terminal)
I found that GDAL and OSSIM agreed well in the forward transform, but while they both had roughly similar sized errors in the inverse they were still fairly distinct results. Reading through the code the reasoning became fairly obvious. They both use iterative solutions for the inverse transform and start with very different initial guesses, and they have fairly weak requirements for covergence. In the case of GDAL the convergence threshold was 1/10th of a pixel. OK, so that leaves me thinking there may be room to improve the inverse functions but accepting that the errors I'm seeing are really all about the iterative convergence.
Next I tried using the RPCs with gdalwarp to actually "rectify" the image. Initially I just left a default elevation of zero for gdalwarp, so I used a command like:
gdalwarp -rpc simulated_rgb.vrt ortho_0.tif
The results seemed plausible but I didn't have a particularly good way to validate it. My test scene is in Indiana, so I used a lat/long to height web service to determine that the region of my scene had a rough elevation of 254m above sea level (we will treat that as the same as the ellipsoid for now). So I tried:
gdalwarp -rpc -to RPC_HEIGHT=254 simulated_rgb.vrt ortho_254.tif
On review of the resulting image, it exactly overlayed the original. I had been hoping the result would have been slightly different based on the elevation. After a moments reflection and a review of the RPC which had many zero values in the terms related to elevation I realized I had't provided sufficient inputs to actually model the perspective transform. I had't given any GCPs at different heights. Also, the example code made it clear that an elevation sensitive RPC could not be generated with less than 80 GCPs (since the equations have 80 unknowns).
My next step was to produce more than 80 GCPs and to ensure that there was more than one elevation represented. It seemed like two layers of 49 GCPs - one at elevation zero, and one at an elevation of 300m should be enough to do the trick, so I updated my GCP generating script to produce an extra set with an ellipsoid that was 300m "bigger". I then ran that through gcps2rpc.
I was pleased to discover that I still got average and maximum errors at the GCPs very similar to the all-zero case. I also confirmed that the GCPs at 300m were actually at a non-trivially different location. So it seemed the RPCs were modelling the perspective view well, at least at the defining positions. I thought I was home free! So I ran a gcpwarp at a height of 0m, 300m and 254m. The results at 0m and 300m looked fine, but somehow the results at 254 meters was radically wrong - greatly inflated! A quick test with gdaltransform gave a sense of what was going wrong:
gdaltransform -rpc simulated_rgb.vrt
0 0 0
-85.036383862713 39.6431112403042 0
0 0 300
-85.0363998100256 39.6431575414329 300
0 0 254
-84.9815303476089 39.2011322215184 254
Essentially this showed that pixel/line location 0,0 was transforming reasonably at the elevations 0 and 300, but at 254 the results were far from the region expected. Clearly the RPCs were not generally modelling the perspective view well. This is the point I am at now. Is it not clear to me what I will need to do to produce stable RPCs for my perspective view from orbit. Perhaps there is a more suitable pattern or number of GCPs to produce. Perhaps the OSSIM RPC generation code needs to be revisited. Perhaps I need to manually construct RPC coefficients that model the perspective view instead of just depending on a least squares best fit to produce reasonable RPCs from a set of GCPs.
Anyways, I hope that this post will provide a bit of background on RPCs to others walking this road, and I hope to have an actually solution to post about in the future.
Sunday, September 8, 2013
Thursday, March 28, 2013
Getting Started with GeoTools
TL;DR: Use Eclipse->Import->Maven to build an eclipse project for GeoTools.
The project I work on at Google (Earth Engine) uses the GeoTools Java library, primarily for coordinate system transformations. We have been using it for a few years for this purpose and it mostly works ok though we have worked around a few quirks. I think it is desirable for our team to be able to contribute at least some small fixes back upstream to GeoTools and so I've been struggling to get to the point where that can be done smoothly.
The next step was internal, where my coworker +Eric Engle, who maintains our internal GeoTool tree, adjusted our use of GeoTools so we build from source and are in a position to patch our local version as we improve things. I had tried to do this once last fall and failed to overcome some challenges related to special requirements of manifest documents for GeoTools. Something to do with finding factories for things like EPSG databases and such.
And now every time I want my geotools clone on github updated I need to issue the following commands locally in my local github clone - basically pulling the changes down to my local repository and then pushing it back up to my own clone:
git fetch upstream
git merge upstream/master
git push origin master
Blech.
I think things would be a bit different if I had write permission on the GeoTools git repository as then I'd be able to make my changes in a branch in it and then send pull requests against that but that would mean lots of branching which I'm also not too fond of.
The default Java on my workstation is a (possibly googlized) OpenJDK 1.7. I get the impression from Jody that GeoTools is targetted at JDK 1.6 and that most folks use the Oracle JDK, not OpenJDK. I was able to able to revert to using OpenJDK 1.6 but was still getting lots of things missing (I wish I had kept better notes on this part). The net result after searching around was that Maven has "profiles" instructing it how to find some classes in different environments that seems to mostly affect the generation of Java docs. After some messing around we found that because my OpenJDK returned "Google Inc." as the vendor something special was needed and Jody was kind enough to incorporate this upstream.
After this change I was able to compile a Maven build of GeoTools successfully, but it took hours and hours for reasons that were never clear and likely somehow Google specific. Today I decided to try this again at home on my fairly stock Ubuntu 10.04 machine. While not exactly fast, it was much faster and I could compile a "mvn install" on an already built tree without testing in a couple minutes. I will note this time I used a downloaded Maven instead of the one provided by Ubuntu, but trying that at work did not seem to help.
The project I work on at Google (Earth Engine) uses the GeoTools Java library, primarily for coordinate system transformations. We have been using it for a few years for this purpose and it mostly works ok though we have worked around a few quirks. I think it is desirable for our team to be able to contribute at least some small fixes back upstream to GeoTools and so I've been struggling to get to the point where that can be done smoothly.
Contribution Agreements
The first challenge was that we needed permission to contribute. Generally this is not a problem as Google is pretty open source friendly, but in order to get the GeoTools corporate contribution agreement signed we needed to clear it with our legal folks. They balked at signing the original agreement which was ... ahem ... hand crafted ... during the GeoTools incubation process. At this point I was prepared to back off on the effort as I didn't want be the guy who raises a stink because his companies legal people were being difficult, as if that should make it the rest of the worlds problems. However, +Jody Garnett was concerned that others might run into similar concerns and pushed through a GeoTools change request adopting a slight variation on the more widely accepted Apache contribution agreement.The next step was internal, where my coworker +Eric Engle, who maintains our internal GeoTool tree, adjusted our use of GeoTools so we build from source and are in a position to patch our local version as we improve things. I had tried to do this once last fall and failed to overcome some challenges related to special requirements of manifest documents for GeoTools. Something to do with finding factories for things like EPSG databases and such.
MapProjection Tolerances
This put things back in my court. I had a small internally reproducable problem with the spherical mercator projection not staying within the error bounds expected by GeoTools at some points on the globe. After some debugging Eclipse I came to the conclusion that:- It is evil that MapProjection.transform() actually uses Java asserts to confirm that point transformation are invertable within a tolerance. This code only triggers when assertions are enabled, and makes it impossible to run code with transformations outside a reasonable realm without triggering this error when assertions are enabled.
- That the default getToleranceForAssertions() result is awfully tight - 1E-5 meters or 1/10th of millimeter!
- For spherical projections it checks orthodromicDistance() which after stepping through the code I determined is very sensitive to rounding error for small distances. Most locations near where I ran into a problem rounded to zero for distance in meters, while the slightly greater distance in radians at the problem point evaluated to something more like one centimeter, exceeding the error threshold.
Well, those are the details, and I hope eventually to prepare a proper GeoTools bug report with a proposed patch. In order to get to that I decided I needed to setup a "regular" GeoTools build that I could run tests against, and use to prepare a patch. The GeoTools manual has great getting started information which I tried to follow.
Git Hate
GeoTools is hosted on github and I have a love/hate relationship with git. Hmm, lets amend that. I have a hate/hate relationship with git. As I somewhat jestingly tell people I haven't yet got a full decade out of subversion yet after learning it in the middle of last decade so I start with a grudge about being forced to learn a new version control technology while the old one continues to work fine from my perspective. But I also find git way more complicated than necessary. Ultimately I determined that I needed to clone geotools on github so I would have a public repository that I could stage changes on and send pull-requests from since I can't issue pull requests for repositories living on my workstations at Google or at home. Then I cloned that repository locally at work and eventually at home. OK, I struggled a bit before concluding that is what was needed but I'm ok with that.
Next I came to the conclusion that there is no easy way from the github interface to refresh my github clone from the official git instance, and I had to setup "upstream" on my workstation clone with the command:
git remote add upstream http://github.com/geotools/geotools.gitAnd now every time I want my geotools clone on github updated I need to issue the following commands locally in my local github clone - basically pulling the changes down to my local repository and then pushing it back up to my own clone:
git fetch upstream
git merge upstream/master
git push origin master
Blech.
I think things would be a bit different if I had write permission on the GeoTools git repository as then I'd be able to make my changes in a branch in it and then send pull requests against that but that would mean lots of branching which I'm also not too fond of.
Maven
GeoTools uses Maven for building. I installed the Maven Ubuntu package and seemed good to go. All I should need to do is "mvn install" in the GeoTools root directory to build.The default Java on my workstation is a (possibly googlized) OpenJDK 1.7. I get the impression from Jody that GeoTools is targetted at JDK 1.6 and that most folks use the Oracle JDK, not OpenJDK. I was able to able to revert to using OpenJDK 1.6 but was still getting lots of things missing (I wish I had kept better notes on this part). The net result after searching around was that Maven has "profiles" instructing it how to find some classes in different environments that seems to mostly affect the generation of Java docs. After some messing around we found that because my OpenJDK returned "Google Inc." as the vendor something special was needed and Jody was kind enough to incorporate this upstream.
After this change I was able to compile a Maven build of GeoTools successfully, but it took hours and hours for reasons that were never clear and likely somehow Google specific. Today I decided to try this again at home on my fairly stock Ubuntu 10.04 machine. While not exactly fast, it was much faster and I could compile a "mvn install" on an already built tree without testing in a couple minutes. I will note this time I used a downloaded Maven instead of the one provided by Ubuntu, but trying that at work did not seem to help.
Eclipse
Next, I wanted to get to the point where I had an Eclipse project I could use to debug and improve GeoTools with. I tried to following the Eclipse Quick Start docs and learned a few things.- Stock Eclipse 3.8 does not seem to include any of the Maven support even though there is apparently some in Eclipse 3.7. I did find downloading and installing the latest Eclipse (Juno - 4.something) did have good Maven support.
- The described process produces a project where you can use GeoTools as jars, and even debug in, but not modify the GeoTools source. This is ok for a use-only GeoTools user but not for someone hoping to fix things or even make temporary modifications during debugging.
I stumbled around quite a bit trying to get past this. What I discovered is that the Eclipse File->Import->Maven->Existing Maven Projects path can be used to setup a project for one or more GeoTools modules. Select the root directory of the GeoTools source tree and Eclipse will walk the pom.xml files (Maven "makefiles) to figure out the dependencies and where the source is. You will need to select which ones to include.
I ran into complaints about lacking Eclipse understanding of git-commit-id-plugin and something about jjtree. I was able to work around this by not including the "cql" module which needed the jjtree stuff, and by hacking out the entire git-commit-id-plugin from the GeoTools root pom.xml file. After that I also restricted myself to selecting the main "library" and I think "ogc" modules.
The result was an Eclipse project in my workspace for each GeoTools module, including the referencing module I was interested in. I was able to go into the src/test tree and run the junit tests fairly easily and all succeeded except three calls to *.class.desiredAssertionStatus() which I commented out. They did not seem terribly relevant.
Conclusion
Well, I can't say that I have accomplished my first meaningful contribution to GeoTools yet, but I have at least gotten to the point where I can build and work on things. Despite some great getting start docs (mostly by Jody I think) and direct support from Jody in IRC (#geotools on freenode) it was still a painful process. This is partly due to my relative lack of familiarity with Maven, Github, and even Eclipse outside of the very narrow way I use it at Google. But it also seems that GeoTools is a big project and there are lots of traps one can fall into with versions of the various tools (Eclipse, Maven, Java runtime, etc). But hopefully people will see that if I can hack it, you likely can too.
Thursday, March 15, 2012
PROJ 4.8.0 Released
Early this week I announced the release of PROJ 4.8.0. It turns out the last release (4.7.0) was in 2009 so it has been quite some time coming.
PROJ.4 is a C library implementation transformations between a wide variety of geospatial coordinate systems. It features include support for many projection methods (ie. Transverse Mercator) as well as support for converting between geodetic datums. As a library it is used by many of the C/C++ open source geospatial applications, and a few proprietary ones including GDAL/OGR, PostGIS, GRASS, MapServer and QGIS (some of my favorites). It also includes commandline utilities mostly used for testing and a dictionary of coordinate system definitions based on the EPSG coordinate system from OGP.
The announcement lists the major new features, but there is a subset I would like to call out here because I think they are quite significant.
I would also be negligent if I failed to call out the important funding support of the US Army Corp of Engineers in these developments via Howard Butler and Liblas / PDAL.
PROJ.4 is a C library implementation transformations between a wide variety of geospatial coordinate systems. It features include support for many projection methods (ie. Transverse Mercator) as well as support for converting between geodetic datums. As a library it is used by many of the C/C++ open source geospatial applications, and a few proprietary ones including GDAL/OGR, PostGIS, GRASS, MapServer and QGIS (some of my favorites). It also includes commandline utilities mostly used for testing and a dictionary of coordinate system definitions based on the EPSG coordinate system from OGP.
The announcement lists the major new features, but there is a subset I would like to call out here because I think they are quite significant.
- A great deal of work went into making the library threadsafe. The main previous hole in threadsafe operation was the use of the global pj_errno for errors. This has been worked around by introduction of an "execution context" object (projCtx) which can be passed in when creating projection objects (projPJ). The execution context also now provides a mechanism to manage error and debug message output. See use of projCtx in proj_api.h for an idea of how this would be used. The execution context also means that more than one part of a process can use PROJ.4 with less risk of stepping on each other (imagine MapServer and GDAL/OGR for instance).
- Preliminary support has been added for vertical datums using the +geoidgrids= directive. The geoidgrids are assumed to provide a transformation from geoid elevations to ellipsoidal elevations and they are discussed a bit in the VerticalDatums topic on the web site.
- I have introduced yet another supported horizontal grid shift file format which I call CTable2. It is really the same as the venerable binary format produced by nad2bin, but with the header regularized and the byte order fixed so that the format is binary compatible across all platforms. This means it is no longer to worry about having distinct copies of the grid shift files for environments like 32bit vs. 64bit linux, or windows vs. linux. I had originally intended to just migrate to using NTv2 but it turns out to take twice as much space because of the presence of rarely used error estimate bands.
- The addition of the +axis= switch to control axis orientation. This makes it possible to handle coordinate systems where "X" is north and "Y" is east as well as various other combinations.
I would also be negligent if I failed to call out the important funding support of the US Army Corp of Engineers in these developments via Howard Butler and Liblas / PDAL.
Saturday, November 5, 2011
FOSS4G Videos
Andrew Ross has announced that all the FOSS4G 2011 and State of the Map videos are now up and available at the FOSSLC web site. Andrew provided the following urls to scan these two categories which gives an extended blurb about each video:
http://www.fosslc.org/drupal/category/event/foss4g2011
http://www.fosslc.org/drupal/category/event/sotm
However I prefer the search mechanism which gives lots more talks per page with just the title.
http://www.fosslc.org/drupal/search/node/foss4g2011
At the risk of being self obsessed I'll highlight two of my talks. I can't seem to find the other two.
I had also really wanted to attend Martin Davis' talk on Spatial Processing Using JEQL but somehow missed it. This weekend I'll review it. There were some other great keynotes and presentations. As you can imagine Paul Ramsey's keynote on Open Source Business Models was informative, amusing and engaging.
The videos were captured by Andrew's team of volunteers, and deployed on the FOSSLC site which includes lots of great videos of talks from other conferences. This has really be a passion of Andrew's for some time and at least in the OSGeo community I don't think the value of it has been fully recognized and exploited. I hope that others will look through the presentation and blog, tweet, +1 and like some that they find the most useful and interesting. Virtual word of mouth is one way to improve the leverage from the work of capturing the videos and the effort to produce the original presentations.
I'd also like to thank Andrew and his helpers Teresa Baldwin, Scott Mitchell, Christine Richard, Alex Trofast, Thanha Ha, "Steve", Paul Deshamps, Assefa Yewondwossen, Nathalie West, Mike Adair and "Ben", many of whom I count as personal friends from the OSGeo Ottawa crowd.
http://www.fosslc.org/drupal/category/event/foss4g2011
http://www.fosslc.org/drupal/category/event/sotm
However I prefer the search mechanism which gives lots more talks per page with just the title.
http://www.fosslc.org/drupal/search/node/foss4g2011
At the risk of being self obsessed I'll highlight two of my talks. I can't seem to find the other two.
I had also really wanted to attend Martin Davis' talk on Spatial Processing Using JEQL but somehow missed it. This weekend I'll review it. There were some other great keynotes and presentations. As you can imagine Paul Ramsey's keynote on Open Source Business Models was informative, amusing and engaging.
The videos were captured by Andrew's team of volunteers, and deployed on the FOSSLC site which includes lots of great videos of talks from other conferences. This has really be a passion of Andrew's for some time and at least in the OSGeo community I don't think the value of it has been fully recognized and exploited. I hope that others will look through the presentation and blog, tweet, +1 and like some that they find the most useful and interesting. Virtual word of mouth is one way to improve the leverage from the work of capturing the videos and the effort to produce the original presentations.
I'd also like to thank Andrew and his helpers Teresa Baldwin, Scott Mitchell, Christine Richard, Alex Trofast, Thanha Ha, "Steve", Paul Deshamps, Assefa Yewondwossen, Nathalie West, Mike Adair and "Ben", many of whom I count as personal friends from the OSGeo Ottawa crowd.
Monday, June 20, 2011
Joining Google
Today I accepted a job with Google as a GIS Data Engineer. I will be based in Mountain View California at head office, and involved in various sorts of geodata processing though I don't really know the details of my responsibilities yet.
I have received occasional email solicitations from Google recruiters in the past, but hadn't really taken them too seriously. I assumed it was a "wide net" search for job candidates and I wasn't looking for regular employment anyways. I was happy enough in my role of independent consultant. It gave me great flexibility and a quite satisfactory income.
Various things conspired to make me contemplate my options. A friend, Phil Vachon, moved to New York City for a fantastic new job. His description of the compensation package made me realize opportunities might exist that I was missing out on. Also, with separation from my wife last year the need to stay in Canada was reduced. So when a Google recruiter contacted me again this spring I took a moment to look over the job description.
The description was for the role of GIS Data Engineer it was a dead on match for my skill set. So, I thought I would at least investigate a bit, and responded. I didn't hear back for many weeks so I assumed I had been winnowed out already or perhaps the position had already been filled. But a few weeks ago I was contacted, and invited to participate in a phone screening interview. That went well, and among other things my knowledge of GIS file formats, and coordinate systems stood me in good stead. So I was invited to California for in person interviews.
That was a grueling day. Five interviews including lots of "write a program to do this on the black board" sorts of questions. Those interviewing me, mostly engineers in the geo group, seemed to know little of my background and I came away feeling somewhat out of place, and like they were really just looking for an entry level engineer. Luckily Michael Weiss-Malik, spent some time at lunch talking about what his group does and gave me a sense of where I would fit in. This made me more comfortable.
Last week I received and offer and it was quite generous. It certainly put my annual consulting income to shame. Also contact with a couple friends inside Google gave me a sense that there were those advocating on my behalf who did know more about my strengths.
I agonized over the weekend and went back and forth quite a bit on the whole prospect. While the financial offer was very good, I didn't particularly need that. And returning to "regular employment" was not inline with my hopes to travel widely, doing my consulting from a variety of exotic locales as long as they had internet access. But, by salting away money it would make it much easier to pursue such a lifestyle at some point in the future.
I also needed to consult extensively with my family. It was very difficult to leave my kids behind, though it helps to know that as teens they are quite independent, I am also confident that their mother would be right there taking good care of them. My kids, and other family members were all very supportive.
My other big concern was giving up my role in the open source geospatial community. While nothing in the job description prevents me from still being active in the GDAL project, and OSGeo, it is clear the job will consume much of my energy and time. I can't expect to play the role I have played in the past, particularly of enabling use of GDAL in commercial contexts by virtual of being available as a paid resource to support users.
I don't have the whole answer to this, but it is certainly my intention to remain active in GDAL/OGR and in OSGeo (on the board and other committees). One big plus of working at Google is the concept of 20% time. I haven't gotten all the details, but it is roughly allowing me to use 20% of my work time for self directed projects. My hope is to use much of my 20% time to work on GDAL/OGR.
Google does make use of GDAL/OGR for some internal data processing and in products like Google Earth Professional. My original hope had been that my job would at least partly be in support of GDAL and possibly other open source technologies within Google. While things are still a bit vague that does not seem to be immediately the case though I'm optimistic such opportunities might arise in the future. But I think this usage does mean that work on GDAL is a reasonable thing to spend 20% time on.
I also made it clear that I still planned to participate in OSGeo events like the FOSS4G conference. I'm pleased to confirm that I will be attending FOSS4G 2011 in Denver as a Google representative and I am confident this will be possible in future years as well.
The coming months will involve many changes for me, and I certainly have a great deal to learn to make myself effective as an employee of Google. But I am optimistic that this will be a job and work place that still allows me to participate in, and contribute to, the community that I love so much. I also think this will be a great opportunity for me to grow. Writing file translators for 20 years can in some ways become a rut!
I have received occasional email solicitations from Google recruiters in the past, but hadn't really taken them too seriously. I assumed it was a "wide net" search for job candidates and I wasn't looking for regular employment anyways. I was happy enough in my role of independent consultant. It gave me great flexibility and a quite satisfactory income.
Various things conspired to make me contemplate my options. A friend, Phil Vachon, moved to New York City for a fantastic new job. His description of the compensation package made me realize opportunities might exist that I was missing out on. Also, with separation from my wife last year the need to stay in Canada was reduced. So when a Google recruiter contacted me again this spring I took a moment to look over the job description.
The description was for the role of GIS Data Engineer it was a dead on match for my skill set. So, I thought I would at least investigate a bit, and responded. I didn't hear back for many weeks so I assumed I had been winnowed out already or perhaps the position had already been filled. But a few weeks ago I was contacted, and invited to participate in a phone screening interview. That went well, and among other things my knowledge of GIS file formats, and coordinate systems stood me in good stead. So I was invited to California for in person interviews.
That was a grueling day. Five interviews including lots of "write a program to do this on the black board" sorts of questions. Those interviewing me, mostly engineers in the geo group, seemed to know little of my background and I came away feeling somewhat out of place, and like they were really just looking for an entry level engineer. Luckily Michael Weiss-Malik, spent some time at lunch talking about what his group does and gave me a sense of where I would fit in. This made me more comfortable.
Last week I received and offer and it was quite generous. It certainly put my annual consulting income to shame. Also contact with a couple friends inside Google gave me a sense that there were those advocating on my behalf who did know more about my strengths.
I agonized over the weekend and went back and forth quite a bit on the whole prospect. While the financial offer was very good, I didn't particularly need that. And returning to "regular employment" was not inline with my hopes to travel widely, doing my consulting from a variety of exotic locales as long as they had internet access. But, by salting away money it would make it much easier to pursue such a lifestyle at some point in the future.
I also needed to consult extensively with my family. It was very difficult to leave my kids behind, though it helps to know that as teens they are quite independent, I am also confident that their mother would be right there taking good care of them. My kids, and other family members were all very supportive.
My other big concern was giving up my role in the open source geospatial community. While nothing in the job description prevents me from still being active in the GDAL project, and OSGeo, it is clear the job will consume much of my energy and time. I can't expect to play the role I have played in the past, particularly of enabling use of GDAL in commercial contexts by virtual of being available as a paid resource to support users.
I don't have the whole answer to this, but it is certainly my intention to remain active in GDAL/OGR and in OSGeo (on the board and other committees). One big plus of working at Google is the concept of 20% time. I haven't gotten all the details, but it is roughly allowing me to use 20% of my work time for self directed projects. My hope is to use much of my 20% time to work on GDAL/OGR.
Google does make use of GDAL/OGR for some internal data processing and in products like Google Earth Professional. My original hope had been that my job would at least partly be in support of GDAL and possibly other open source technologies within Google. While things are still a bit vague that does not seem to be immediately the case though I'm optimistic such opportunities might arise in the future. But I think this usage does mean that work on GDAL is a reasonable thing to spend 20% time on.
I also made it clear that I still planned to participate in OSGeo events like the FOSS4G conference. I'm pleased to confirm that I will be attending FOSS4G 2011 in Denver as a Google representative and I am confident this will be possible in future years as well.
The coming months will involve many changes for me, and I certainly have a great deal to learn to make myself effective as an employee of Google. But I am optimistic that this will be a job and work place that still allows me to participate in, and contribute to, the community that I love so much. I also think this will be a great opportunity for me to grow. Writing file translators for 20 years can in some ways become a rut!
Friday, February 18, 2011
MapServer TIFF Overview Performance
Last week Thomas Bonfort, MapServer polymath contacted me with some surprising performance results he was getting while comparing two ways of handling image overviews with MapServer.
In one case he had a single GeoTIFF file with overviews built by GDAL's gdaladdo utility in a separate .ovr file. Only one MapServer LAYER was used to refer to this image. In the second case he used a distinct TIFF file for each power-of-two overview level, and a distinct LAYER with minscale/maxscale settings for each file. So in the first case it was up to GDAL to select the optimum overview level out of a merged file, while in the second case MapServer decided which overview to use via the layer scale selection logic.
He was seeing performance using the MapServer multi-layer approach being about twice what it was with letting GDAL select the right overview. He prepared a graph and scripts with everything I needed to reproduce his results.
I was able to reproduce similar results at my end, and set to work trying to analyse what was going on. The scripts used apache-bench to run MapServer in it's usual cgi-bin incarnation. This is a typical use case, and in aggregate shows the performance. However, my usual technique for investigating performance bottlenecks is very low-tech. I do a long run demonstrating the issue in gdb, and hit cntl-c frequently, and use "where" to examine what is going on. If the bottleneck is dramatic enough this is usually informative. But I could not apply this to the very very short running cgi-bin processes.
So I set out to establish a similar work load in a single long running MapScript application. In doing so I discovered a few things. First, that there was no way to easily load the OWSRequest parameters from an url except through the QUERY_STRING environment variable. So I extended mapscript/swiginc/owsrequest.i to have a loadParamsFromURL() method. My test script then looked like:
The second thing I learned is that Thomas' recent work to directly use libjpeg and libpng for output in MapServer had not honoured the msIO_ IO redirection mechanism needed for the above. I fixed that too.
This gave me a proces that would run for a while and that I could debug with gdb. A sampling of "what is going on" showed that much of the time was being spent loading TIFF tags from directories - particularly the tile offset and tile size tag values.
The base file used is 130000 x 150000 pixels, and is internally broken up into nearly 900000 256x256 tiles. Any map request would only use a few tiles but the array of pointers to tiles and their sizes for just the base image amounted to approximately 14MB. So in order to get about 100K of imagery out of the files we were also reading at least 14MB of tile index/sizes.
The GDAL overview case was worse than the MapServer overview case because when opening the file GDAL scans through all the overview directories to identify what overviews are available. This means we have to load the tile offset/size values for all the overviews regardless of whether we will use them later. When the offset/size values are read to scan a directory they are subsequently discarded when the next directory is read. So in cases where we need to come back to a particular overview level we still have to reload the offsets and sizes.
For several years I have had some concerns about the efficiency of files with large tile offset/size arrays, and with the cost of jumping back and forth between different overviews with GDAL. This case highlighted the issue. It also suggested an obvious optimization - to avoid loading the tile offset and size values until we actually need them. If we access an image directory (ie overview level) to get general information about it such as the size, but we don't actually need imagery we could completely skip reading the offsets and sizes.
So I set about implementing this in libtiff. It is helpful being a core maintainer of libtiff as well as GDAL and Mapserver. :-) The change was a bit messy, and also it seemed a bit risky due to how libtiff handled directory tags. So I treat the change as as an experimental (default off) build time option controlled by the new libtiff configure option --enable-defer-strile-load.
With the change in place, I now get comparable results using the two different overview strategies. Yipee! Pleased with this result I have commited the changes in libtiff CVS head, and pulled them downstream into the "built in" copy of libtiff in GDAL.
However, it occurs to me that there is still an unfortunate amount of work being done to load the tile offset/size vectors when doing full resolution image accesses for MapServer. A format that computed the tile offset and size instead of storing all the offsets and sizes explicitly might do noticeable better for this use case. In fact, Erdas Imagine (with a .ige spill file) is such a format. Perhaps tonight I can compare to that and contemplate a specialized form of uncompressed TIFF files which doesn't need to store and load the tile offset/sizes.
I would like to thank Thomas Bonfort for identifying this issue, and providing such an easy to use set of scripts to reproduce and demonstrate the problem. I would also like to thank MapGears and USACE for supporting my time on this and other MapServer raster problems.
In one case he had a single GeoTIFF file with overviews built by GDAL's gdaladdo utility in a separate .ovr file. Only one MapServer LAYER was used to refer to this image. In the second case he used a distinct TIFF file for each power-of-two overview level, and a distinct LAYER with minscale/maxscale settings for each file. So in the first case it was up to GDAL to select the optimum overview level out of a merged file, while in the second case MapServer decided which overview to use via the layer scale selection logic.
He was seeing performance using the MapServer multi-layer approach being about twice what it was with letting GDAL select the right overview. He prepared a graph and scripts with everything I needed to reproduce his results.
I was able to reproduce similar results at my end, and set to work trying to analyse what was going on. The scripts used apache-bench to run MapServer in it's usual cgi-bin incarnation. This is a typical use case, and in aggregate shows the performance. However, my usual technique for investigating performance bottlenecks is very low-tech. I do a long run demonstrating the issue in gdb, and hit cntl-c frequently, and use "where" to examine what is going on. If the bottleneck is dramatic enough this is usually informative. But I could not apply this to the very very short running cgi-bin processes.
So I set out to establish a similar work load in a single long running MapScript application. In doing so I discovered a few things. First, that there was no way to easily load the OWSRequest parameters from an url except through the QUERY_STRING environment variable. So I extended mapscript/swiginc/owsrequest.i to have a loadParamsFromURL() method. My test script then looked like:
import mapscript
map = mapscript.mapObj( 'cnes.map')
mapscript.msIO_installStdoutToBuffer()
for i in range(1000):
req = mapscript.OWSRequest()
req.loadParamsFromURL( 'LAYERS=truemarble-gdal&FORMAT=image/jpeg&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&STYLES=&EXCEPTIONS=application/vnd.ogc.se_inimage&SRS=EPSG%3A900913&BBOX=1663269.7343875,1203424.5723063,1673053.6740063,1213208.511925&WIDTH=256&HEIGHT=256')
map.OWSDispatch( req )
The second thing I learned is that Thomas' recent work to directly use libjpeg and libpng for output in MapServer had not honoured the msIO_ IO redirection mechanism needed for the above. I fixed that too.
This gave me a proces that would run for a while and that I could debug with gdb. A sampling of "what is going on" showed that much of the time was being spent loading TIFF tags from directories - particularly the tile offset and tile size tag values.
The base file used is 130000 x 150000 pixels, and is internally broken up into nearly 900000 256x256 tiles. Any map request would only use a few tiles but the array of pointers to tiles and their sizes for just the base image amounted to approximately 14MB. So in order to get about 100K of imagery out of the files we were also reading at least 14MB of tile index/sizes.
The GDAL overview case was worse than the MapServer overview case because when opening the file GDAL scans through all the overview directories to identify what overviews are available. This means we have to load the tile offset/size values for all the overviews regardless of whether we will use them later. When the offset/size values are read to scan a directory they are subsequently discarded when the next directory is read. So in cases where we need to come back to a particular overview level we still have to reload the offsets and sizes.
For several years I have had some concerns about the efficiency of files with large tile offset/size arrays, and with the cost of jumping back and forth between different overviews with GDAL. This case highlighted the issue. It also suggested an obvious optimization - to avoid loading the tile offset and size values until we actually need them. If we access an image directory (ie overview level) to get general information about it such as the size, but we don't actually need imagery we could completely skip reading the offsets and sizes.
So I set about implementing this in libtiff. It is helpful being a core maintainer of libtiff as well as GDAL and Mapserver. :-) The change was a bit messy, and also it seemed a bit risky due to how libtiff handled directory tags. So I treat the change as as an experimental (default off) build time option controlled by the new libtiff configure option --enable-defer-strile-load.
With the change in place, I now get comparable results using the two different overview strategies. Yipee! Pleased with this result I have commited the changes in libtiff CVS head, and pulled them downstream into the "built in" copy of libtiff in GDAL.
However, it occurs to me that there is still an unfortunate amount of work being done to load the tile offset/size vectors when doing full resolution image accesses for MapServer. A format that computed the tile offset and size instead of storing all the offsets and sizes explicitly might do noticeable better for this use case. In fact, Erdas Imagine (with a .ige spill file) is such a format. Perhaps tonight I can compare to that and contemplate a specialized form of uncompressed TIFF files which doesn't need to store and load the tile offset/sizes.
I would like to thank Thomas Bonfort for identifying this issue, and providing such an easy to use set of scripts to reproduce and demonstrate the problem. I would also like to thank MapGears and USACE for supporting my time on this and other MapServer raster problems.
Saturday, December 18, 2010
RPATH, RUNPATH and LD_LIBRARY_PATH
Long time since my last post - my bad.
Recently I have been working on a renewed FWTools for linux - hopefully I'll follow up on that later. Much like in the past FWTools is intended to work across a variety of linux systems (ie. Ubuntu, Redhat, Suse, Gentoo) and to not interfere with the native packing system nor to require special permissions to install.
To that end we try to distribute as many of the required shared libraries as possible, and we use the LD_LIBRARY_PATH environment variable to ensure FWTools finds and uses these shared libraries in preference to any system provided versions which might be incompatible (though often aren't).
In my local testing I was getting frustrated as my shared libraries were not being used. Instead the system versions of things like libsqlite3.so.4 were getting picked up from /usr/lib64 despite my LD_LIBRARY_PATH pointing to a copy distributed with FWTools. I check this with the ldd command which shows what sublibraries an executable or shared library need and which is being picked according to the current policy.
I carefully read the man page for ld.so, which is responsible for loading shared libraries at runtime, and it indicated:
The necessary shared libraries needed by the program are searched for
in the following order
o Using the environment variable LD_LIBRARY_PATH
(LD_AOUT_LIBRARY_PATH for a.out programs). Except if the exe‐
cutable is a setuid/setgid binary, in which case it is ignored.
...
Well, my programs are not setuid/setgid nor are they a.out (very old style of linux binary), so what gives? Due to poor time synchronization on my systems and use of nfs, I noticed that some shared library dates were in the future. I thought perhaps some funky security policy prevented them from being used, but that theory did not pan out. I tried doing it on a local file system incase there were some rule about not trusting shared libraries from remote file systems in some cases. Still no go.
However, googling around the net I found that my local ld.so man page was not very comprehensive. Another ld.so man page gave me a hint, referring to DT_RPATH and DT_RUNPATH attributes being used before LD_LIBRARY_PATH. I was vaguely aware of shared libraries having an "rpath" attribute which could be set at link time and might have some sort of influence on searching. Distribution packaging folks have sometimes bugged me about it, especially back in the old days before libtool was used for GDAL. But I had never really understood it well and it just seemed like unnecessary stuff.
Further digging confirmed that ELF (the modern linux object file format) executables and shared libraries may have an RPATH attribute built right into the file. This is essentially a directory that should be searched first for shared libraries that this executable or shared library depends on. RUNPATH is similar, but is apparently evaluated after LD_LIBRARY_PATH. I got my critical information from this blog post by Trevor King.
My old favorite, the ldd command, does not show the rpath settings though it does honour them. I also use nm to dump the symbols in object files and shared libraries, but it doesn't address special attributes like rpath. Digging around I found the readelf command which can report all sorts of information about a shared library or executable, including the RPATH, and RUNPATH. It also helpfully lists the shared libraries that are needed, something that often gets obscured in ldd by the chaining of dependencies.
The rpath attribute is generally set by libtool when linking shared libraries. GDAL includes a build option to build without libtool, and I contemplated using that which would help alot. But it wouldn't help for other libraries in FWTools that were not directly pulled in by GDAL. My goal with the modern FWTools is to - as much as possible - use standard CentOS/RedHat rpm binaries as the basis for the product. So I don't want to have to rebuild everything from scratch. What I really wanted was a way to suppress use of the rpath, switch it to runpath or remove it entirely. Low and behold, I am not the first with this requirement. There is a command called chrpath that can do any of these "in place" to a binary (see chrpath(1) man page).
This command isn't installed by default on my system, but it was packaged and it works very nicely. I'll preprocess the shared libraries and binaries to strip out rpath so my LD_LIBRARY_PATH magic works smoothly. I had even considered using it to replace the rpath at install time with the path inside the FWTools distribution for libraries. This, unfortunately, does not appear to be possible as it seems chrpath can only replace rpath's with rpaths of the same or a shorter length. In any event, FWTools still depends on cover scripts to set other environment variables so doing this wouldn't have made it practical to avoid my environment setup magic.
I owe a lot to some other blog and mailing list posts turned up by google. This blog entry is partly an attempt to seed the net with helpful answers for others that run into similar problems. I will say it is very annoying that a company decided to name itself "rPath" - it made searching quite a bit harder!
Recently I have been working on a renewed FWTools for linux - hopefully I'll follow up on that later. Much like in the past FWTools is intended to work across a variety of linux systems (ie. Ubuntu, Redhat, Suse, Gentoo) and to not interfere with the native packing system nor to require special permissions to install.
To that end we try to distribute as many of the required shared libraries as possible, and we use the LD_LIBRARY_PATH environment variable to ensure FWTools finds and uses these shared libraries in preference to any system provided versions which might be incompatible (though often aren't).
In my local testing I was getting frustrated as my shared libraries were not being used. Instead the system versions of things like libsqlite3.so.4 were getting picked up from /usr/lib64 despite my LD_LIBRARY_PATH pointing to a copy distributed with FWTools. I check this with the ldd command which shows what sublibraries an executable or shared library need and which is being picked according to the current policy.
I carefully read the man page for ld.so, which is responsible for loading shared libraries at runtime, and it indicated:
The necessary shared libraries needed by the program are searched for
in the following order
o Using the environment variable LD_LIBRARY_PATH
(LD_AOUT_LIBRARY_PATH for a.out programs). Except if the exe‐
cutable is a setuid/setgid binary, in which case it is ignored.
...
Well, my programs are not setuid/setgid nor are they a.out (very old style of linux binary), so what gives? Due to poor time synchronization on my systems and use of nfs, I noticed that some shared library dates were in the future. I thought perhaps some funky security policy prevented them from being used, but that theory did not pan out. I tried doing it on a local file system incase there were some rule about not trusting shared libraries from remote file systems in some cases. Still no go.
However, googling around the net I found that my local ld.so man page was not very comprehensive. Another ld.so man page gave me a hint, referring to DT_RPATH and DT_RUNPATH attributes being used before LD_LIBRARY_PATH. I was vaguely aware of shared libraries having an "rpath" attribute which could be set at link time and might have some sort of influence on searching. Distribution packaging folks have sometimes bugged me about it, especially back in the old days before libtool was used for GDAL. But I had never really understood it well and it just seemed like unnecessary stuff.
Further digging confirmed that ELF (the modern linux object file format) executables and shared libraries may have an RPATH attribute built right into the file. This is essentially a directory that should be searched first for shared libraries that this executable or shared library depends on. RUNPATH is similar, but is apparently evaluated after LD_LIBRARY_PATH. I got my critical information from this blog post by Trevor King.
My old favorite, the ldd command, does not show the rpath settings though it does honour them. I also use nm to dump the symbols in object files and shared libraries, but it doesn't address special attributes like rpath. Digging around I found the readelf command which can report all sorts of information about a shared library or executable, including the RPATH, and RUNPATH. It also helpfully lists the shared libraries that are needed, something that often gets obscured in ldd by the chaining of dependencies.
The rpath attribute is generally set by libtool when linking shared libraries. GDAL includes a build option to build without libtool, and I contemplated using that which would help alot. But it wouldn't help for other libraries in FWTools that were not directly pulled in by GDAL. My goal with the modern FWTools is to - as much as possible - use standard CentOS/RedHat rpm binaries as the basis for the product. So I don't want to have to rebuild everything from scratch. What I really wanted was a way to suppress use of the rpath, switch it to runpath or remove it entirely. Low and behold, I am not the first with this requirement. There is a command called chrpath that can do any of these "in place" to a binary (see chrpath(1) man page).
chrpath -d usr/lib64/libgdal.so
This command isn't installed by default on my system, but it was packaged and it works very nicely. I'll preprocess the shared libraries and binaries to strip out rpath so my LD_LIBRARY_PATH magic works smoothly. I had even considered using it to replace the rpath at install time with the path inside the FWTools distribution for libraries. This, unfortunately, does not appear to be possible as it seems chrpath can only replace rpath's with rpaths of the same or a shorter length. In any event, FWTools still depends on cover scripts to set other environment variables so doing this wouldn't have made it practical to avoid my environment setup magic.
I owe a lot to some other blog and mailing list posts turned up by google. This blog entry is partly an attempt to seed the net with helpful answers for others that run into similar problems. I will say it is very annoying that a company decided to name itself "rPath" - it made searching quite a bit harder!
Sunday, July 18, 2010
OGR DXF Upgrade
Well, I just finished up a week of work implementing PCIDSK Vector write and update support. My next task is an upgrade to the OGR DXF driver requested by Stadt Uster to better meet their production requirements.
They need the ability to control layer naming, production of dash patterns, and producing objects as block references.
Currently the DXF writer is dependent on everything except the entities section coming directly from a template file. This has meant it was not practical to create layers, line styles and block references on the fly. This is going to have to change now, though we will continue to use the template header extensively.
The first approach considered was just to require the user to develop a template header with all the layer names, line styles and block references predefined. This might have been adequate for Uster who have specific needs and once a template was prepared they could generally just reuse it for additional products. But it would have made the new capabilities of very little utility to other users of GDAL.
So the planned approach has two parts. First we scan the template header for layer definitions, and block definitions (and possibly line style definitions). Then as we go through the entities if we find these layers or blocks referenced we just use them directly.
However, if we find layers referenced from the objects being written to the DXF file (based on the "Layer" attribute) we will automatically create a new layer in the header, matching the configuration of the default layer ("0"). This means we can automatically create layers that have identity even if they are otherwise indistinguishable.
For block references we use a roughly similar approach. We prescan for block definitions, but then we extend them with any entities written to an OGR layer named "Blocks". Normally all DXF entities are exposed through an OGR layer called "entities" though when writing we accept any layer name, except now for Blocks which is special. Then when we write we allow these blocks to be referenced based on a BlockName attribute.
Corresponding behavior will also be available when reading. If desired (a config option turned on) we will expose block definitions as a Blocks layer, and the actual block references in the entities layer will just be a point feature with insertion information. The default behavior will remain what it does not - which is to inline copies of the block geometries for each block reference as this is the only approach that will be handled gracefully by most applications or writers.
For line styles it is not clear that it is helpful to predefine them (though I'll examine that during implementation). The only aspect we are interested in preserving and producing is dash-dot patterns. But while writing I will create new line styles for each such pattern encountered, and then write them out to the header.
Currently the DXF writer copies the template header to the output file, then writes entities, then appends the trailer template. In the future I will need to write the entities to a temporary file, storing up block, layer and line style definitions in memory. Then on closing the dataset I will need to compose the full header, append the entities and trailer.
I dislike this pattern. It introduces a need for temporary files which can lead to surprising disk use requirements. Also, if something goes wrong the temp files may not get cleaned up properly. We also lose any hope of streaming operation. However, to achieve what we want to with the DXF driver it seems unavoidable.
Hopefully I'll have the DXF changes in trunk within a couple weeks. If there are folks interested in DXF generation, keep an eye on SVN for updates. Beta testers appreciated!
They need the ability to control layer naming, production of dash patterns, and producing objects as block references.
Currently the DXF writer is dependent on everything except the entities section coming directly from a template file. This has meant it was not practical to create layers, line styles and block references on the fly. This is going to have to change now, though we will continue to use the template header extensively.
The first approach considered was just to require the user to develop a template header with all the layer names, line styles and block references predefined. This might have been adequate for Uster who have specific needs and once a template was prepared they could generally just reuse it for additional products. But it would have made the new capabilities of very little utility to other users of GDAL.
So the planned approach has two parts. First we scan the template header for layer definitions, and block definitions (and possibly line style definitions). Then as we go through the entities if we find these layers or blocks referenced we just use them directly.
However, if we find layers referenced from the objects being written to the DXF file (based on the "Layer" attribute) we will automatically create a new layer in the header, matching the configuration of the default layer ("0"). This means we can automatically create layers that have identity even if they are otherwise indistinguishable.
For block references we use a roughly similar approach. We prescan for block definitions, but then we extend them with any entities written to an OGR layer named "Blocks". Normally all DXF entities are exposed through an OGR layer called "entities" though when writing we accept any layer name, except now for Blocks which is special. Then when we write we allow these blocks to be referenced based on a BlockName attribute.
Corresponding behavior will also be available when reading. If desired (a config option turned on) we will expose block definitions as a Blocks layer, and the actual block references in the entities layer will just be a point feature with insertion information. The default behavior will remain what it does not - which is to inline copies of the block geometries for each block reference as this is the only approach that will be handled gracefully by most applications or writers.
For line styles it is not clear that it is helpful to predefine them (though I'll examine that during implementation). The only aspect we are interested in preserving and producing is dash-dot patterns. But while writing I will create new line styles for each such pattern encountered, and then write them out to the header.
Currently the DXF writer copies the template header to the output file, then writes entities, then appends the trailer template. In the future I will need to write the entities to a temporary file, storing up block, layer and line style definitions in memory. Then on closing the dataset I will need to compose the full header, append the entities and trailer.
I dislike this pattern. It introduces a need for temporary files which can lead to surprising disk use requirements. Also, if something goes wrong the temp files may not get cleaned up properly. We also lose any hope of streaming operation. However, to achieve what we want to with the DXF driver it seems unavoidable.
Hopefully I'll have the DXF changes in trunk within a couple weeks. If there are folks interested in DXF generation, keep an eye on SVN for updates. Beta testers appreciated!
Thursday, June 17, 2010
Travelling
I have spent this week in beautiful Chicoutimi / Saguenay this week at the Rendez-vous OSGeo Québec. This is a two day conference put on by the local chapter here in Quebec in association with another local geomatics conference.
It has been some time since I have been to a significant event, unfortunately having missed FOSS4G in Sydney, so it was good to be involved again. As long ago as the Ottawa OSGIS conference in 2004 I recall Daniel Morissette, Steve Lime, and me hanging out in Dave McIllhagga's penthouse and talking about visiting Chicoutimi. Perhaps because of that Daniel, who was a primary organizer, made an effort to invite Steve and me. Unfortunately Steve was unable to attend in person due to a knee injury, but many other friends were here including Paul Ramsey and Jeff McKenna. As you can imagine, in addition to an excellent program, there were great social activities.
Today, during the other geomatics conference, I am helping to man the OSGeo booth. I unfortunately speak no french, and I was worried about having nothing I could say or provide to french-only speakers. However, I was immensely pleased to discover that the booth was well equipted with the full set of OSGeo promotional materials all translated into french, printed out and very professional looking.
I have at times not been all that excited about translation efforts. I am uni-lingual, and can't really help. I have also had the impression that at least some of the translation efforts have been rather haphazard and might not reflected particularly well on the projects or OSGeo. But here I was very impressed with the complete and professional job that has been done with the promotional materials, and also how useful it is. It really is very helpful for our credibility today, and I can see this would apply in many other venues. I understand the translation is done by the Francophone chapter, and that one major contributor has been Yves Jacolin. I am sure there are many others also involved. So, my hat is off to their efforts.
I should also mention, I finally got to see Paul's fantastic talk "Beyond Nerds Bearing Gifts: The future of the Open Source Economy". I can see why everyone was raving about it after Sydney. One message that I took away from the talk was that as all companies become open source companies to some degree, we need to be cautious about thinking of ourselves primarily as open source companies or consultants. Our excellent integration into the open source ecosystem is one of our selling points, but we need to ensure our primary focus is on serving the needs of our clients effectively and efficiently.
Today I was also invited to attend a technical meeting in Zurich at the end of the month as part of the GEM (Global Earthquate Model) project. I was not aware of this activity at all until Steve Arnold (aka nerdboy on IRC) brought it to my attention recently. Among other things, GEM is interested in using open source components in their IT infrastructure and I will provide some feedback on their plans and efforts so far, particularly on the geospatial aspects.
This week, after discussions with Paul Ramsey, Jeff McKenna, and Venkatesh Raghavan, I am tentatively planning to go to Japan in November for the OSGeo meetings held there. Everyone who has attended these meetings in the past as an invited speaker has raved about what a good experience it is, so I am looking forward to it very much. I believe that Paul Ramsey will also likely be there, and Jeff may be back in Japan at the same time for his OSGeo4W cooperation.
Also, everyone is very much looking forward to Barcelona. I am confident it will be the best FOSS4G ever. I have registered, but still need to book my flight and hotel. I'm thinking of taking a little extra time to see the area. I sometimes find travel tiring and puts me behind in my work. But I think the social and networking opportunities are too good to miss. These meetings reengerize me in ways that last much longer than the jet lag.
It has been some time since I have been to a significant event, unfortunately having missed FOSS4G in Sydney, so it was good to be involved again. As long ago as the Ottawa OSGIS conference in 2004 I recall Daniel Morissette, Steve Lime, and me hanging out in Dave McIllhagga's penthouse and talking about visiting Chicoutimi. Perhaps because of that Daniel, who was a primary organizer, made an effort to invite Steve and me. Unfortunately Steve was unable to attend in person due to a knee injury, but many other friends were here including Paul Ramsey and Jeff McKenna. As you can imagine, in addition to an excellent program, there were great social activities.
Today, during the other geomatics conference, I am helping to man the OSGeo booth. I unfortunately speak no french, and I was worried about having nothing I could say or provide to french-only speakers. However, I was immensely pleased to discover that the booth was well equipted with the full set of OSGeo promotional materials all translated into french, printed out and very professional looking.
I have at times not been all that excited about translation efforts. I am uni-lingual, and can't really help. I have also had the impression that at least some of the translation efforts have been rather haphazard and might not reflected particularly well on the projects or OSGeo. But here I was very impressed with the complete and professional job that has been done with the promotional materials, and also how useful it is. It really is very helpful for our credibility today, and I can see this would apply in many other venues. I understand the translation is done by the Francophone chapter, and that one major contributor has been Yves Jacolin. I am sure there are many others also involved. So, my hat is off to their efforts.
I should also mention, I finally got to see Paul's fantastic talk "Beyond Nerds Bearing Gifts: The future of the Open Source Economy". I can see why everyone was raving about it after Sydney. One message that I took away from the talk was that as all companies become open source companies to some degree, we need to be cautious about thinking of ourselves primarily as open source companies or consultants. Our excellent integration into the open source ecosystem is one of our selling points, but we need to ensure our primary focus is on serving the needs of our clients effectively and efficiently.
Today I was also invited to attend a technical meeting in Zurich at the end of the month as part of the GEM (Global Earthquate Model) project. I was not aware of this activity at all until Steve Arnold (aka nerdboy on IRC) brought it to my attention recently. Among other things, GEM is interested in using open source components in their IT infrastructure and I will provide some feedback on their plans and efforts so far, particularly on the geospatial aspects.
This week, after discussions with Paul Ramsey, Jeff McKenna, and Venkatesh Raghavan, I am tentatively planning to go to Japan in November for the OSGeo meetings held there. Everyone who has attended these meetings in the past as an invited speaker has raved about what a good experience it is, so I am looking forward to it very much. I believe that Paul Ramsey will also likely be there, and Jeff may be back in Japan at the same time for his OSGeo4W cooperation.
Also, everyone is very much looking forward to Barcelona. I am confident it will be the best FOSS4G ever. I have registered, but still need to book my flight and hotel. I'm thinking of taking a little extra time to see the area. I sometimes find travel tiring and puts me behind in my work. But I think the social and networking opportunities are too good to miss. These meetings reengerize me in ways that last much longer than the jet lag.
Sunday, March 14, 2010
Slaying the Datum Shift Dragon
In the last few weeks, I believe I have made substantial improvements to the way datum shift information is derived from the EPSG dictionary for use in GDAL/OGR, PROJ.4 and related packages.
The EPSG database model supports having multiple datum shift options for a particular datum, such as "Potsdam Rauenberg 1950 DHDN" to get to WGS84. This reflects the reality that datum conversion is often an approximation for which there may be multiple reasonable solutions. Often different datum shift parameters are appropriate depending on the sub-region of the datum being used.
Many other coordinate system dictionary models and GIS systems do *not* provide for a multiplicity of datum shift options. For instance, the OGC Well Known Text representation of a coordinate system only allows for one TOWGS84[] clause indicating the preferred mapping to WGS84. This is also true of many software packages, including those based on PROJ.4 and it's epsg based dictionary.
In the past, the code I had implemented for translating the epsg database to a dictionary format assumed that datum shift parameters should only be carried into the dictionary if there was one, and only one shift available in the EPSG database. I took this approach because it seemed very dangerous to arbitrarily select one of several possible datum shifts without any intrinsic knowledge about which was most appropriate. I depending on users seeing that no datum shift was available in the default translation and taking this to mean that they had to do some research to establish what was most appropriate for them.
Predictably, the real result was massive confusion and complaints. At one time if there was no datum shift available, PROJ.4 would just do a transformation based on the change of ellipsoid which was often a very poor choice. So in PROJ 4.6.0 I altered the code not attempt any datum shift transformation if either or both of the source and destination coordinate system lacked datum shift information.
While this got rid of one family of errors, it also triggered lots of additional frustration and confusion. Part of this was just because there was a change of behavior. But it was also pushing people a bit harder to determine the appropriate datum shift and they did not like having to do this.
So, at last, I have taken the plunge and reworked the scripts used to translate the EPSG dictionary so that they attempt to pick one shift if there are several available. I follow a few heuristics in this effort.
1) I discard any datum shifts marked as deprecated.
2) I examine the supersession table to identify any datum shifts that have been superceeded by newer forms and ignore the superceeded forms.
3) I try to pick the datum shift with the larges "area of use" region under the assumption that it will likely be the broad use shift rather than a shift only applicable in a small area.
4) I examine the datum_shift_pref.csv file to see if there is a user supplied preferred datum shift to use. If so, I use that.
The result of all this is a datum_shift.csv file which includes all the datum shifts, with one of them marked as preferred. That preferred version goes into the gcs.csv file for use with the associated geographic coordinate system.
This seems to be working reasonable well, and I did a big pass through open tickets in GDAL, libgeotiff and PROJ.4 to find outstanding datum shift issues. I believe that the bulk are now resolved.
Currently GDAL, and the PROJ.4 dictionary are just using the preferred datum shift. But the intention of keeping the datum_shift.csv file with all the possible shifts for any given GCS is that savvy applications could let the user choose. Also, my hope is that an each mechanism will be added in the future so that users can alter the preferred setting in the datum_shift.csv file, and then have the gcs.csv regenerated. That waits to be done.
I'm also fixing a few other issues in the coordinate system realm including support for axis orientation (ie. South Orientated Transverse Mercator) and fixing a few other translations. I'd like to thank INGRES who have supported this work as part of an effort to bring top notch coordinate system support into the INGRES geospatial project. I'd also like to thank Jan Hartmann, and Mikael Rittri who assisted with suggested approaches, and verification of the results.
The EPSG database model supports having multiple datum shift options for a particular datum, such as "Potsdam Rauenberg 1950 DHDN" to get to WGS84. This reflects the reality that datum conversion is often an approximation for which there may be multiple reasonable solutions. Often different datum shift parameters are appropriate depending on the sub-region of the datum being used.
Many other coordinate system dictionary models and GIS systems do *not* provide for a multiplicity of datum shift options. For instance, the OGC Well Known Text representation of a coordinate system only allows for one TOWGS84[] clause indicating the preferred mapping to WGS84. This is also true of many software packages, including those based on PROJ.4 and it's epsg based dictionary.
In the past, the code I had implemented for translating the epsg database to a dictionary format assumed that datum shift parameters should only be carried into the dictionary if there was one, and only one shift available in the EPSG database. I took this approach because it seemed very dangerous to arbitrarily select one of several possible datum shifts without any intrinsic knowledge about which was most appropriate. I depending on users seeing that no datum shift was available in the default translation and taking this to mean that they had to do some research to establish what was most appropriate for them.
Predictably, the real result was massive confusion and complaints. At one time if there was no datum shift available, PROJ.4 would just do a transformation based on the change of ellipsoid which was often a very poor choice. So in PROJ 4.6.0 I altered the code not attempt any datum shift transformation if either or both of the source and destination coordinate system lacked datum shift information.
While this got rid of one family of errors, it also triggered lots of additional frustration and confusion. Part of this was just because there was a change of behavior. But it was also pushing people a bit harder to determine the appropriate datum shift and they did not like having to do this.
So, at last, I have taken the plunge and reworked the scripts used to translate the EPSG dictionary so that they attempt to pick one shift if there are several available. I follow a few heuristics in this effort.
1) I discard any datum shifts marked as deprecated.
2) I examine the supersession table to identify any datum shifts that have been superceeded by newer forms and ignore the superceeded forms.
3) I try to pick the datum shift with the larges "area of use" region under the assumption that it will likely be the broad use shift rather than a shift only applicable in a small area.
4) I examine the datum_shift_pref.csv file to see if there is a user supplied preferred datum shift to use. If so, I use that.
The result of all this is a datum_shift.csv file which includes all the datum shifts, with one of them marked as preferred. That preferred version goes into the gcs.csv file for use with the associated geographic coordinate system.
This seems to be working reasonable well, and I did a big pass through open tickets in GDAL, libgeotiff and PROJ.4 to find outstanding datum shift issues. I believe that the bulk are now resolved.
Currently GDAL, and the PROJ.4 dictionary are just using the preferred datum shift. But the intention of keeping the datum_shift.csv file with all the possible shifts for any given GCS is that savvy applications could let the user choose. Also, my hope is that an each mechanism will be added in the future so that users can alter the preferred setting in the datum_shift.csv file, and then have the gcs.csv regenerated. That waits to be done.
I'm also fixing a few other issues in the coordinate system realm including support for axis orientation (ie. South Orientated Transverse Mercator) and fixing a few other translations. I'd like to thank INGRES who have supported this work as part of an effort to bring top notch coordinate system support into the INGRES geospatial project. I'd also like to thank Jan Hartmann, and Mikael Rittri who assisted with suggested approaches, and verification of the results.
Wednesday, February 10, 2010
GDAL 1.7.0 is dead, long live GDAL 1.7.1
Due to a serious bug in GDAL 1.7.0 resulting in all Erdas Imagine files produced by the release being unreadable to any non-GDAL applications, including Erdas Imagine, the project has decided to retract the 1.7.0 release. In it's place today we have issued a GDAL/OGR 1.7.1 release. I sincerely hope that use of 1.7.0 in the wild will not end up producing unreadable Erdas Imagine files that haunt us for years.
The introduction of the bug was due to a very subtle issue, and it didn't affect GDAL's ability to read the files. So our extensive regression tests for this file format did not detect any problems. But the lesson would appear to be that we still need mechanisms to do cross testing with other packages as often as possible. We are looking into ways of more vigorous and organized cross testing at major releases.
The introduction of the bug was due to a very subtle issue, and it didn't affect GDAL's ability to read the files. So our extensive regression tests for this file format did not detect any problems. But the lesson would appear to be that we still need mechanisms to do cross testing with other packages as often as possible. We are looking into ways of more vigorous and organized cross testing at major releases.
Tuesday, February 9, 2010
World Mapping
Over much of the last decade I have fought with world mapping issues in MapServer. This falls into three broad areas - problems with requests crossing the dateline, problems with requests including the poles, and problems with requests extending beyond the domain of validity of a projection. Typically these are not a big issue when working on local mapping projects, but do become issues when doing world level views, or when exploring areas far from the domain of validity of some of the data.
Several years ago I made a modest effort to collect some hints and tips, and knowledge about the issues into a topic in the old MapServer wiki. With the planned decommissioning of the old mapserver web site at UMN (and in light of recent work in this area), I have migrated the old wiki topic to the MapServer Trac Wiki as OldWorldMappingIssues, and started an updated wiki topic on the issues.
The problems I have been investigating in recent months on behalf of a client, are related to providing WMS service in web mercator projection with MapServer, but cascading to another WMS providing data in EPSG:4326 (geographic). The initial problem report was with WMS requests in the Alaska/Siberia area that were oddly pixelated. The problem was that the request region from -21000000 to -17000000 easting was getting reprojected to geographic as -179 to 179 in longitude. This was occurring because value further east than the dateline (roughly -20000000 meters) were getting returned as positive longitudes since PROJ.4 returns longitude values in the range -180 to 180. This transformed BBOX was sent to the cascaded WMS, but with a width and height that reflected the approximate size of the original request. The result was that the map pulled from the cascaded WMS had very poor horizontal resolution (since 360 degrees of longitude were being reflected in relatively few pixels instead of the roughly 10 degrees that should have been fetched).
The solution I developed for this as part of ticket #3179 was to utilize MapServer's special "dateline line/polygon wrapping logic" in msProjectRect(). This logic basically makes it so that small polygons crossing the dateline will be kept local even if the longitudes go outside the -180 to 180 default range. The special dateline wrapping logic was implemented several years ago but was only used when reprojecting shapes. This worked well for this specific problem, but I avoided incorporating it into MapServer 5.6 since I was concerned it might be quite a disruptive change to drop in near release time.
Last week my client contacted me with a related problem. The default view of the world is roughly -21000000,-11000000,21000000,11000000 in web mercator. This is slightly beyond the real extent of the world, so by default PROJ.4 maps -2100000 easting to around +179 longitude even though I might prefer to have thought of it as -181 longitude. The simplest transformation of the web mercator BBOX to lat/long gives an area of roughly -179 to 179 longitude, and -90 to 90 latitude. This is a reasonably close approximation, but loses a bit of data at the edges. This is an issue that has come up in several forms in the past and is sometimes is helped by using of the LOAD_WHOLE_IMAGE processing option.
But my client was running into a different problem. They were seeing WMS requests with a BBOX of -179.79,-89.00,323.96,85.07. The funny item being the longitude of 323. It turns out that this is due to relatively poor handling of world scale polygons in the dateline wrapping logic. It rewraps some points and then realizes that the polygon is really not a small local polygon and it leaves some points unaltered. This can result, as here, in polygons more than 360 degrees wide. This causes some sort of flipping anomaly which I did not investigate.
After some analysis and hair pulling I came to the conclusion that it was essentially impossible to make the polygon dateline wrapping robust in the face of mercator to geographic transformations that had already wrapped things into the -180 to 180 domain. Even if I backed out the approach entirely, I'd still end up with substantial over-requesting in some cases. For this particular problem I did institute a bit of logic in msProjectRect() that said a geographic rect that is wider than 360 degrees should be remapped as minx=-180 and maxx=180. But I was able to construct other requests that would still behave badly.
What I *really* wanted was for PROJ.4 to reproject mercator values of less than -20000000 easting as less than -180 longitude, and more than 20000000 easting as more than 180 longitude. In desperation I finally delved into PROJ.4 to better understand how this wrapping occurs and if their might be some way of avoiding it. I was flabbergasted to discover PROJ.4 was doing the -180/180 wrapping in pj_inv() and that there was a +over switch that could be used to turn this wrapping behavior off!
For my client, I am now advising they use +over in their web mercator projection definition, and this helps avoid many dateline related reprojection problems. I've also added some documentation on this in the PROJ.4 General Parameters documentation. (BTW, a shout out to Hamish who has improved this doc quite a bit!)
I'm planning to incorporate my reproject rectangle as a polygon logic in MapServer trunk and I hope people will keep it in mind as a potential source of other problems. But I think wider knowledge of, and use of +over may be the big win of this whole effort! I also hope people will contribute to the World Mapping Issues topic in the MapServer wiki.
Several years ago I made a modest effort to collect some hints and tips, and knowledge about the issues into a topic in the old MapServer wiki. With the planned decommissioning of the old mapserver web site at UMN (and in light of recent work in this area), I have migrated the old wiki topic to the MapServer Trac Wiki as OldWorldMappingIssues, and started an updated wiki topic on the issues.
The problems I have been investigating in recent months on behalf of a client, are related to providing WMS service in web mercator projection with MapServer, but cascading to another WMS providing data in EPSG:4326 (geographic). The initial problem report was with WMS requests in the Alaska/Siberia area that were oddly pixelated. The problem was that the request region from -21000000 to -17000000 easting was getting reprojected to geographic as -179 to 179 in longitude. This was occurring because value further east than the dateline (roughly -20000000 meters) were getting returned as positive longitudes since PROJ.4 returns longitude values in the range -180 to 180. This transformed BBOX was sent to the cascaded WMS, but with a width and height that reflected the approximate size of the original request. The result was that the map pulled from the cascaded WMS had very poor horizontal resolution (since 360 degrees of longitude were being reflected in relatively few pixels instead of the roughly 10 degrees that should have been fetched).
The solution I developed for this as part of ticket #3179 was to utilize MapServer's special "dateline line/polygon wrapping logic" in msProjectRect(). This logic basically makes it so that small polygons crossing the dateline will be kept local even if the longitudes go outside the -180 to 180 default range. The special dateline wrapping logic was implemented several years ago but was only used when reprojecting shapes. This worked well for this specific problem, but I avoided incorporating it into MapServer 5.6 since I was concerned it might be quite a disruptive change to drop in near release time.
Last week my client contacted me with a related problem. The default view of the world is roughly -21000000,-11000000,21000000,11000000 in web mercator. This is slightly beyond the real extent of the world, so by default PROJ.4 maps -2100000 easting to around +179 longitude even though I might prefer to have thought of it as -181 longitude. The simplest transformation of the web mercator BBOX to lat/long gives an area of roughly -179 to 179 longitude, and -90 to 90 latitude. This is a reasonably close approximation, but loses a bit of data at the edges. This is an issue that has come up in several forms in the past and is sometimes is helped by using of the LOAD_WHOLE_IMAGE processing option.
But my client was running into a different problem. They were seeing WMS requests with a BBOX of -179.79,-89.00,323.96,85.07. The funny item being the longitude of 323. It turns out that this is due to relatively poor handling of world scale polygons in the dateline wrapping logic. It rewraps some points and then realizes that the polygon is really not a small local polygon and it leaves some points unaltered. This can result, as here, in polygons more than 360 degrees wide. This causes some sort of flipping anomaly which I did not investigate.
After some analysis and hair pulling I came to the conclusion that it was essentially impossible to make the polygon dateline wrapping robust in the face of mercator to geographic transformations that had already wrapped things into the -180 to 180 domain. Even if I backed out the approach entirely, I'd still end up with substantial over-requesting in some cases. For this particular problem I did institute a bit of logic in msProjectRect() that said a geographic rect that is wider than 360 degrees should be remapped as minx=-180 and maxx=180. But I was able to construct other requests that would still behave badly.
What I *really* wanted was for PROJ.4 to reproject mercator values of less than -20000000 easting as less than -180 longitude, and more than 20000000 easting as more than 180 longitude. In desperation I finally delved into PROJ.4 to better understand how this wrapping occurs and if their might be some way of avoiding it. I was flabbergasted to discover PROJ.4 was doing the -180/180 wrapping in pj_inv() and that there was a +over switch that could be used to turn this wrapping behavior off!
For my client, I am now advising they use +over in their web mercator projection definition, and this helps avoid many dateline related reprojection problems. I've also added some documentation on this in the PROJ.4 General Parameters documentation. (BTW, a shout out to Hamish who has improved this doc quite a bit!)
PROJECTION
"+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378137 +b=6378137
+units=m +over +no_defs"
END
I'm planning to incorporate my reproject rectangle as a polygon logic in MapServer trunk and I hope people will keep it in mind as a potential source of other problems. But I think wider knowledge of, and use of +over may be the big win of this whole effort! I also hope people will contribute to the World Mapping Issues topic in the MapServer wiki.
Friday, January 15, 2010
Visiting China This Month
In a diversion from my usual programming related geo-geeking, I'm heading to China on January 20th to do some in person, on the ground examination of geography. I will be visiting a friend there, and returning on the 29th. I will be visiting Shanghai, and Fuzhou (on the coast roughly half way between Shanghai and Hong Kong).
Today I received my visa from the Chinese embassy. This is the first time I've ever had to get a visa for international travel, and I have a new respect for those who have to do it regularly to attend events like FOSS4G. It is nearly a five hour activity to travel to Ottawa and visit the embassy, which I had to do to drop off my request, and to pick up the visa. Also, they were going to refuse to accept my request until I could provide a clear address I would be at. I had to rush out to an internet cafe to book a hotel for the first night just so I would have an address to write on the form. I barely got back before the office closed (9am to 1pm!), which would have cost me another day of travel. Normally I prefer to play things by ear after I arrive.
I'm flying via Newark despite my preference to avoid going through the USA given the security hassles that entails. But the price difference over flying Air Canada through Vancouver was just too significant.
I'm not sure what my connectivity will be, so I'm hoping to load up my laptop with some work I can do offline if needed. Likely I will focus on write/update access for PCIDSK vectors in the PCIDSK SDK.
Today I received my visa from the Chinese embassy. This is the first time I've ever had to get a visa for international travel, and I have a new respect for those who have to do it regularly to attend events like FOSS4G. It is nearly a five hour activity to travel to Ottawa and visit the embassy, which I had to do to drop off my request, and to pick up the visa. Also, they were going to refuse to accept my request until I could provide a clear address I would be at. I had to rush out to an internet cafe to book a hotel for the first night just so I would have an address to write on the form. I barely got back before the office closed (9am to 1pm!), which would have cost me another day of travel. Normally I prefer to play things by ear after I arrive.
I'm flying via Newark despite my preference to avoid going through the USA given the security hassles that entails. But the price difference over flying Air Canada through Vancouver was just too significant.
I'm not sure what my connectivity will be, so I'm hoping to load up my laptop with some work I can do offline if needed. Likely I will focus on write/update access for PCIDSK vectors in the PCIDSK SDK.
Monday, January 4, 2010
hsv_merge.py
Recently Trent Hare indicated an interest in an extension to the gdaldem program to produce a color relief image mixed with a hillshade.
The gdaldem utility is a utility program for producing various products from DEMs (elevation models). It was originally written by Matthew Perry, with contributions by Even Rouault, Howard Butler, and Chris Yesson to incorporate it into GDAL proper. It can produce color relief maps, hillshades, slope, aspect and various roughness values.
But what Trent wanted was a product with both a hillshade and a color relief image merged. The color relief shows overall elevation with color, and the hillshade gives a sense of local structure. He has been doing this with ImageMagick to produce the product from the hillshade and relief produced by gdaldem, but this resulted in the georeferencing being lost. It is accomplished by transforming the relief image to HSV (hue, saturation and value) color space, and then replacing the "value" portion with the hillshade greyscale image. The value is really the intensity or brightness. This is a common technique and is often used to merge high resolution greyscale imagery with lower resolution RGB image to produce a result with meaningful color and high resolution.
I reviewed the gdaldem.cpp code with the thought of adding a new mode, but it turns out the program is fairly complicated and a new hillshade/relief mode would likely be hard to maintain properly. So instead I set out to create a utility that could do the rgb/greyscale merge in hsv space as a standalone utility. When possible I prefer to do non-core utilities as python instead of C++. It feels lighter and easier for folks to munge for other purposes.
In Python we generally try to do image processing using the Python numpy package, so that is what I set out to do. For some things numpy lets us treat arrays as if they were scalars, and operate on the whole array with a single statement. So the "scalar" function to convert hsv to rgb looks like:
I was able to change this into a numpy based function that looked fairly similar in parts, but somewhat different in other areas:
The first step just breaks a single hsv array into three component arrays and is just bookkeeping:
The next few statements look very much like the original and are essentially simple math that happens to be done on whole arrays instead of single values. The only odd part is using "astype(int)" to convert to int instead of the int() builtin function that was used for scalars.
But following this there are many if/then else statements which had looked like:
Unfortunately we can't really do something quite like an if statement with numpy that will be evaluated on a per-pixel (array entry) basis. The old scalar logic assumed that if a condition was met, it could do something for that case, and return. But for arrays we are going to have to do things quite differently. The above is essentially a case statement, so we instead use the numpy choose() function like this:
The choose function goes through each entry in the "i" array, and uses it's value as an index into the set of possible return values (v,q,p,p,t,v in the first case). So if i[0] is 0, then the first item v[0] is assigned to r[0]. If i[0] was 4 then t[0] would be assigned to r[0]. Thus the six if statements are handled as a choose() with six possible return results, and the choose is evaluated for each array entry (pixel).
This works reasonably well for this very structured case. We then finish by merging the red, green and blue values into an 3 band array:
One thing missing in my rendition of this function is the saturation == 0 case:
It is very difficult to have the pixels for which the saturation is zero handled differently than all the rest since we can't actually return control to the caller for only some pixels (array entries). In this particular situation it seems that this was just an optimization and if s is zero the result will be v,v,v anyways so I just omitted this statement.
But this helps highlight the difficulty of taking different logic paths for different entries in the array using numpy. The rgb_to_hsv() logic was even more involved. The original was:
This became:
The key mechanism here is the computation of true/value masks (maxc_is_r, maxc_is_g, maxc_is_b) which are used with choose() statements to overlay different results on the pixels that match particular conditions. It may not be immediately evident, but we are essentially evaluating all the equations for all pixels even though the values are only going to be used for some of them. Quite wasteful.
The final result, hsv_merge.py, works nicely but I am once again left frustrated with the unnatural handling of logic paths in numpy. It took me over three hours to construct this small script. Most of that time was spent poring over the numpy reference trying to figure out how to do things that are dead easy in scalar python.
I love Python, but to be honest I still can't say the same for numpy.
The gdaldem utility is a utility program for producing various products from DEMs (elevation models). It was originally written by Matthew Perry, with contributions by Even Rouault, Howard Butler, and Chris Yesson to incorporate it into GDAL proper. It can produce color relief maps, hillshades, slope, aspect and various roughness values.
But what Trent wanted was a product with both a hillshade and a color relief image merged. The color relief shows overall elevation with color, and the hillshade gives a sense of local structure. He has been doing this with ImageMagick to produce the product from the hillshade and relief produced by gdaldem, but this resulted in the georeferencing being lost. It is accomplished by transforming the relief image to HSV (hue, saturation and value) color space, and then replacing the "value" portion with the hillshade greyscale image. The value is really the intensity or brightness. This is a common technique and is often used to merge high resolution greyscale imagery with lower resolution RGB image to produce a result with meaningful color and high resolution.
I reviewed the gdaldem.cpp code with the thought of adding a new mode, but it turns out the program is fairly complicated and a new hillshade/relief mode would likely be hard to maintain properly. So instead I set out to create a utility that could do the rgb/greyscale merge in hsv space as a standalone utility. When possible I prefer to do non-core utilities as python instead of C++. It feels lighter and easier for folks to munge for other purposes.
In Python we generally try to do image processing using the Python numpy package, so that is what I set out to do. For some things numpy lets us treat arrays as if they were scalars, and operate on the whole array with a single statement. So the "scalar" function to convert hsv to rgb looks like:
def hsv_to_rgb(h, s, v):
if s == 0.0: return v, v, v
i = int(h*6.0) # XXX assume int() truncates!
f = (h*6.0) - i
p = v*(1.0 - s)
q = v*(1.0 - s*f)
t = v*(1.0 - s*(1.0-f))
if i%6 == 0: return v, t, p
if i == 1: return q, v, p
if i == 2: return p, v, t
if i == 3: return p, q, v
if i == 4: return t, p, v
if i == 5: return v, p, q
# Cannot get here
I was able to change this into a numpy based function that looked fairly similar in parts, but somewhat different in other areas:
def hsv_to_rgb( hsv ):
h = hsv[0]
s = hsv[1]
v = hsv[2]
i = (h*6.0).astype(int)
f = (h*6.0) - i
p = v*(1.0 - s)
q = v*(1.0 - s*f)
t = v*(1.0 - s*(1.0-f))
r = i.choose( v, q, p, p, t, v )
g = i.choose( t, v, v, q, p, p )
b = i.choose( p, p, t, v, v, q )
rgb = numpy.asarray([r,g,b]).astype(numpy.byte)
return rgb
The first step just breaks a single hsv array into three component arrays and is just bookkeeping:
h = hsv[0]
s = hsv[1]
v = hsv[2]
The next few statements look very much like the original and are essentially simple math that happens to be done on whole arrays instead of single values. The only odd part is using "astype(int)" to convert to int instead of the int() builtin function that was used for scalars.
i = (h*6.0).astype(int)
f = (h*6.0) - i
p = v*(1.0 - s)
q = v*(1.0 - s*f)
t = v*(1.0 - s*(1.0-f))
But following this there are many if/then else statements which had looked like:
if i%6 == 0: return v, t, p
if i == 1: return q, v, p
if i == 2: return p, v, t
if i == 3: return p, q, v
if i == 4: return t, p, v
if i == 5: return v, p, q
Unfortunately we can't really do something quite like an if statement with numpy that will be evaluated on a per-pixel (array entry) basis. The old scalar logic assumed that if a condition was met, it could do something for that case, and return. But for arrays we are going to have to do things quite differently. The above is essentially a case statement, so we instead use the numpy choose() function like this:
r = i.choose( v, q, p, p, t, v )
g = i.choose( t, v, v, q, p, p )
b = i.choose( p, p, t, v, v, q )
The choose function goes through each entry in the "i" array, and uses it's value as an index into the set of possible return values (v,q,p,p,t,v in the first case). So if i[0] is 0, then the first item v[0] is assigned to r[0]. If i[0] was 4 then t[0] would be assigned to r[0]. Thus the six if statements are handled as a choose() with six possible return results, and the choose is evaluated for each array entry (pixel).
This works reasonably well for this very structured case. We then finish by merging the red, green and blue values into an 3 band array:
rgb = numpy.asarray([r,g,b]).astype(numpy.byte)
One thing missing in my rendition of this function is the saturation == 0 case:
if s == 0.0: return v, v, v
It is very difficult to have the pixels for which the saturation is zero handled differently than all the rest since we can't actually return control to the caller for only some pixels (array entries). In this particular situation it seems that this was just an optimization and if s is zero the result will be v,v,v anyways so I just omitted this statement.
But this helps highlight the difficulty of taking different logic paths for different entries in the array using numpy. The rgb_to_hsv() logic was even more involved. The original was:
def rgb_to_hsv(r, g, b):
maxc = max(r, g, b)
minc = min(r, g, b)
v = maxc
if minc == maxc: return 0.0, 0.0, v
s = (maxc-minc) / maxc
rc = (maxc-r) / (maxc-minc)
gc = (maxc-g) / (maxc-minc)
bc = (maxc-b) / (maxc-minc)
if r == maxc: h = bc-gc
elif g == maxc: h = 2.0+rc-bc
else: h = 4.0+gc-rc
h = (h/6.0) % 1.0
return h, s, v
This became:
def rgb_to_hsv( rgb ):
r = rgb[0]
g = rgb[1]
b = rgb[2]
maxc = numpy.maximum(r,numpy.maximum(g,b))
minc = numpy.minimum(r,numpy.minimum(g,b))
v = maxc
minc_eq_maxc = numpy.equal(minc,maxc)
# compute the difference, but reset zeros to ones to avoid divide by zeros later.
ones = numpy.ones((r.shape[0],r.shape[1]))
maxc_minus_minc = numpy.choose( minc_eq_maxc, (ones, maxc-minc) )
s = (maxc-minc) / numpy.maximum(ones,maxc)
rc = (maxc-r) / maxc_minus_minc
gc = (maxc-g) / maxc_minus_minc
bc = (maxc-b) / maxc_minus_minc
maxc_is_r = numpy.equal(maxc,r)
maxc_is_g = numpy.equal(maxc,g)
maxc_is_b = numpy.equal(maxc,b)
h = numpy.zeros((r.shape[0],r.shape[1]))
h = numpy.choose( maxc_is_b, (h,4.0+gc-rc) )
h = numpy.choose( maxc_is_g, (h,2.0+rc-bc) )
h = numpy.choose( maxc_is_r, (h,bc-gc) )
h = numpy.mod(h/6.0,1.0)
hsv = numpy.asarray([h,s,v])
return hsv
The key mechanism here is the computation of true/value masks (maxc_is_r, maxc_is_g, maxc_is_b) which are used with choose() statements to overlay different results on the pixels that match particular conditions. It may not be immediately evident, but we are essentially evaluating all the equations for all pixels even though the values are only going to be used for some of them. Quite wasteful.
The final result, hsv_merge.py, works nicely but I am once again left frustrated with the unnatural handling of logic paths in numpy. It took me over three hours to construct this small script. Most of that time was spent poring over the numpy reference trying to figure out how to do things that are dead easy in scalar python.
I love Python, but to be honest I still can't say the same for numpy.
Tuesday, December 15, 2009
Death by Complexity
This post will be incomprehensible.
Over the last three days I have spent 4.5 hours working on GDAL Ticket 3276 related to a failure to use external overviews with JPEG2000 compressed NITF files - in particular for NITF files containing more than one jpeg2000 image.
This bug was particularly hairy because it comes at the intersection of several things that are messy/complex in GDAL:
It turns out that .aux.xml metadata was supported for NITF subdatasets, and it was possible to build overviews for nitf subdatasets, and it was possible to substitute external tiff overviews for jpeg2000 data streams in an NITF file. But it was not possible to substitute external tiff overviews for jpeg2000 data stream in an nitf file with multiple images (subdatasets).
It took me a long time to figure out what aspects were broken, and how the various components were supposed to work even though I had implemented most of them. The problem is that many of these capabilities are rarely used, don't fit the standard GDAL model well, and were individually quite complex to implement. The complexity as the various aspects come together is compounded.
Each of the capabilities was added for fairly good reasons - mostly in order to provide a seamless, and performant user experience for GDAL users. But in order to provide this consistent external set of behaviors we are having to build more and more complexity into parts of GDAL - to the point where I am not sure it is sustainable.
Interestingly, most of this complexity has grown without input from the broader GDAL community. The PAM mechanism predated the modern Project Steering Committee and it's RFC process. I added the changes for PAM on subdatasets, and some of the specific NITF driver capabilities without discussion with the PSC on the assumption that they are either bug fixes, or are sufficiently driver-local that the PSC does not need to be involved. Possibly if these changes had needed to be justified in public, push back on the complexity might have prevented some.
The specific problem itself was fixed, as is documented in the ticket, with changeset 18312 holding the core fix. However, this fix is (IMHO) just adding additional fragility to the existing house of cards.
I don't really have a solution to the growing complexity, but perhaps thinking about, and starting to open up the issues a bit is a first step to containing the danger. There is certainly a cautionary tale or two in here.
BTW, GDAL 1.7.0Beta 1 is now released - testing and bug reports are welcome!
Over the last three days I have spent 4.5 hours working on GDAL Ticket 3276 related to a failure to use external overviews with JPEG2000 compressed NITF files - in particular for NITF files containing more than one jpeg2000 image.
This bug was particularly hairy because it comes at the intersection of several things that are messy/complex in GDAL:
- JPEG2000 in NITF is implemented by creating JP2KAK driver dataset wrapping the jpeg2000 image data within the NITF file, and then using it's bands indirectly as bands for the NITF dataset. These band objects still mostly think of the jpeg2000 dataset as "their" dataset but for some purposes we really might wish they knew about the NITF dataset that appears to own them.
- NITF files can contain more than one image. Such multi-image files are treated as containing subdatasets, one per image. For the most part these subdatasets are intended to act as freestanding things, but they are also, to some extent related back to the single file on disk containing the subdatasets.
- Overviews in GDAL are mostly handled through an overview manager object embedded in the GDALDataset base class. However, JPEG2000 images have built-in overviews not handled through the overview manager.
- The PAM (Persistant Auxilary Metadata) mechanism is used via an intermediate GDALPamDataset class to provide a way of storing additional information about datasets that the intrinsic format does not support. This information is stored in an .aux.xml file associated with the main data file.
- In GDAL 1.7 a new capability was added to store PAM information for subdatasets in an .aux.xml file associated with the main data file so that subdatasets would work as much like a regular dataset as possible.
- In GDAL 1.7 support was added for building overviews on subdatsets. Since normally overviews would be stored in a .ovr file with the same basename as the main filename, it was necessary to do something special so that overviews of subdataset would have one .ovr file per subdataset. This was accomplished by keeping the overview file name in the .aux.xml file associated with the subdataset.
It turns out that .aux.xml metadata was supported for NITF subdatasets, and it was possible to build overviews for nitf subdatasets, and it was possible to substitute external tiff overviews for jpeg2000 data streams in an NITF file. But it was not possible to substitute external tiff overviews for jpeg2000 data stream in an nitf file with multiple images (subdatasets).
It took me a long time to figure out what aspects were broken, and how the various components were supposed to work even though I had implemented most of them. The problem is that many of these capabilities are rarely used, don't fit the standard GDAL model well, and were individually quite complex to implement. The complexity as the various aspects come together is compounded.
Each of the capabilities was added for fairly good reasons - mostly in order to provide a seamless, and performant user experience for GDAL users. But in order to provide this consistent external set of behaviors we are having to build more and more complexity into parts of GDAL - to the point where I am not sure it is sustainable.
Interestingly, most of this complexity has grown without input from the broader GDAL community. The PAM mechanism predated the modern Project Steering Committee and it's RFC process. I added the changes for PAM on subdatasets, and some of the specific NITF driver capabilities without discussion with the PSC on the assumption that they are either bug fixes, or are sufficiently driver-local that the PSC does not need to be involved. Possibly if these changes had needed to be justified in public, push back on the complexity might have prevented some.
The specific problem itself was fixed, as is documented in the ticket, with changeset 18312 holding the core fix. However, this fix is (IMHO) just adding additional fragility to the existing house of cards.
I don't really have a solution to the growing complexity, but perhaps thinking about, and starting to open up the issues a bit is a first step to containing the danger. There is certainly a cautionary tale or two in here.
BTW, GDAL 1.7.0Beta 1 is now released - testing and bug reports are welcome!
Saturday, December 5, 2009
OGR DXF Driver
The last couple weeks I have been working on an OGR DXF driver. AutoCAD DXF format is a popular interchange format for CAD data, and to a somewhat lesser extent for goespatial map data often originating from engineering departments. It is a rather ugly format. Even though it is ASCII it is less than fun for humans to scan.
The raw machinery of the format is published by Autodesk, and lots of translators have been written for it in the past. However, I find it very frustrating that the format specifications fail to address the semantics of the format to any meaningful degree. It is assumed, I guess, that the person reading them is already deeply familiar with the AutoCAD data model.
So, for instance, it talks about the BLOCKS section, and the INSERT entity, but it never really explains that by defining a bunch of entities as a block, and then putting them into the drawing it makes it possible to treat a bunch of entities as a group. This has to be deduced. I have an old dxf file produced by FME that uses this as a mechanism to group a set of polylines that form a multipolygon, but I'm not sure how widely this approach is used.
In my driver, I have implemented two mechanisms for supporting blocks. In the case of non-text entities I try to accumulate the geometries I derive from the entities in the block into one OGC geometrycollection. But I also found that blocks in some drawings are used to group sets of text elements. In this case if I only use the geometries, I've discarded the important component - the text - so I collect the text entities as full features. And when I encounter the INSERT entity I push the original set of features into the input stream.
I spent most of a day implementing read support for the DIMENSION entity. This is essentially supposed to be a single object that shows a measurement on a drawing. The red part below is a single dimension element as rendered by QCAD.
Of course, there is no clear analog to this in OGR. So I ended up creating a single feature with a MULTILINESTRING geometry for the leaders, and arrows, and then pushed the label as a separate point feature with a LABEL style string.
I have made some effort to capture the drawing style information so that features can, in theory, be drawn similarly to AutoCAD by OGR applications. There is still really only one way to do this in OGR and that is to provide styling information as OGR Feature Style strings. This is a format that Daniel Morissette, myself and I believe Stephane Villeneuve defined many years ago as based closely on the sorts of styling supported by the mapinfo format(s). A few extensions have been made over the years, and a few OGR drivers utilize it to return styling information including DGN, and MapInfo. Andrey Kiselev has also done some stuff with it, so I presume there is a driver of his that uses it. So far the only two application that I know of using the feature style information are OpenEV and MapServer (via OGR autostyling).
I'm quite conflicted about the OGR feature style specification. I just hate to define "yet another" way of describing feature styling that is not closely related to any proper standard. I can't help but think there ought to be something more standards oriented to use as a basis. Perhaps OGC SLD, or SVG or something like that. But SLD didn't exist when we started, and does not express a lot of what I want. I don't really know that much about SVG - perhaps it would be a good choice. So, lacking a clear plan each time I need styling I do a bit more work on the OGR feature style specification while remaining less than fully committed to it as a long term thing.
Another aspect that was challenging was all the curves. I added the OGRGeometryFactory::approximateArcAngles() method as a generic mechanism to approximate arcs on an ellipsoid or circle as linestrings. The code was adapted from the Oracle driver and similar code exists in the DGN, and NTF drivers too I believe. So, finally this moves into the core.
One of the DXF curve types is a spline. Review the QCAD source code I found they implemented the spline rendering using rationale b-spline curves derived from Chapter 4 of An Introduction to NURBS by David Rogers. They release under the GPL so I can't directly use the code from QCAD without also putting GDAL/OGR under GPL restrictions, so I contacted the original author for permission to publish this code under the MIT/X license, like GDAL. He has indicated he has received my email and is considering the request ... so I wait.
In addition to the specification, I have also found the QCAD (and it's underlying dxflib library) to be a useful reference. I considered using dxflib but it really does not seem to solve any of the hard parts for me, and it would have added a dependency on an external library - complicating building of GDAL/OGR. The other very helpful resource was the v.in.dxf code from GRASS. This is a slightly less sophisticated dxf reader than OGR (IMHO) but it provided a very easy to understand implementation of a DXF reader for GIS data that was especially helpful in writing the proposal for the DXF driver. I also used it as a reference when implementing some of the element translations.
Before closing, I would like to thank Andreas Neumann and the City of Uster who have provided funding for this development. It seems that local government in Switzerland punches above it's weight in the free gis software world (the Kanton of Solothurn is also a big supporter of qgis, and related technologies as I understand). If only one in ten cities in the world made serious use of free gis software and provided enough financial support for one core developer it would have a huge boosting effect.
Anyways, preliminary DXF read support is in GDAL SVN now. Consider trying it out and providing feedbac.
The raw machinery of the format is published by Autodesk, and lots of translators have been written for it in the past. However, I find it very frustrating that the format specifications fail to address the semantics of the format to any meaningful degree. It is assumed, I guess, that the person reading them is already deeply familiar with the AutoCAD data model.
So, for instance, it talks about the BLOCKS section, and the INSERT entity, but it never really explains that by defining a bunch of entities as a block, and then putting them into the drawing it makes it possible to treat a bunch of entities as a group. This has to be deduced. I have an old dxf file produced by FME that uses this as a mechanism to group a set of polylines that form a multipolygon, but I'm not sure how widely this approach is used.
In my driver, I have implemented two mechanisms for supporting blocks. In the case of non-text entities I try to accumulate the geometries I derive from the entities in the block into one OGC geometrycollection. But I also found that blocks in some drawings are used to group sets of text elements. In this case if I only use the geometries, I've discarded the important component - the text - so I collect the text entities as full features. And when I encounter the INSERT entity I push the original set of features into the input stream.
I spent most of a day implementing read support for the DIMENSION entity. This is essentially supposed to be a single object that shows a measurement on a drawing. The red part below is a single dimension element as rendered by QCAD.
Of course, there is no clear analog to this in OGR. So I ended up creating a single feature with a MULTILINESTRING geometry for the leaders, and arrows, and then pushed the label as a separate point feature with a LABEL style string.
I have made some effort to capture the drawing style information so that features can, in theory, be drawn similarly to AutoCAD by OGR applications. There is still really only one way to do this in OGR and that is to provide styling information as OGR Feature Style strings. This is a format that Daniel Morissette, myself and I believe Stephane Villeneuve defined many years ago as based closely on the sorts of styling supported by the mapinfo format(s). A few extensions have been made over the years, and a few OGR drivers utilize it to return styling information including DGN, and MapInfo. Andrey Kiselev has also done some stuff with it, so I presume there is a driver of his that uses it. So far the only two application that I know of using the feature style information are OpenEV and MapServer (via OGR autostyling).
I'm quite conflicted about the OGR feature style specification. I just hate to define "yet another" way of describing feature styling that is not closely related to any proper standard. I can't help but think there ought to be something more standards oriented to use as a basis. Perhaps OGC SLD, or SVG or something like that. But SLD didn't exist when we started, and does not express a lot of what I want. I don't really know that much about SVG - perhaps it would be a good choice. So, lacking a clear plan each time I need styling I do a bit more work on the OGR feature style specification while remaining less than fully committed to it as a long term thing.
Another aspect that was challenging was all the curves. I added the OGRGeometryFactory::approximateArcAngles() method as a generic mechanism to approximate arcs on an ellipsoid or circle as linestrings. The code was adapted from the Oracle driver and similar code exists in the DGN, and NTF drivers too I believe. So, finally this moves into the core.
One of the DXF curve types is a spline. Review the QCAD source code I found they implemented the spline rendering using rationale b-spline curves derived from Chapter 4 of An Introduction to NURBS by David Rogers. They release under the GPL so I can't directly use the code from QCAD without also putting GDAL/OGR under GPL restrictions, so I contacted the original author for permission to publish this code under the MIT/X license, like GDAL. He has indicated he has received my email and is considering the request ... so I wait.
In addition to the specification, I have also found the QCAD (and it's underlying dxflib library) to be a useful reference. I considered using dxflib but it really does not seem to solve any of the hard parts for me, and it would have added a dependency on an external library - complicating building of GDAL/OGR. The other very helpful resource was the v.in.dxf code from GRASS. This is a slightly less sophisticated dxf reader than OGR (IMHO) but it provided a very easy to understand implementation of a DXF reader for GIS data that was especially helpful in writing the proposal for the DXF driver. I also used it as a reference when implementing some of the element translations.
Before closing, I would like to thank Andreas Neumann and the City of Uster who have provided funding for this development. It seems that local government in Switzerland punches above it's weight in the free gis software world (the Kanton of Solothurn is also a big supporter of qgis, and related technologies as I understand). If only one in ten cities in the world made serious use of free gis software and provided enough financial support for one core developer it would have a huge boosting effect.
Anyways, preliminary DXF read support is in GDAL SVN now. Consider trying it out and providing feedbac.
Subscribe to:
Posts (Atom)