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.
Thursday, June 17, 2010
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!
Subscribe to:
Posts (Atom)