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.
Wednesday, February 10, 2010
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)