Global Mapper v25.0

DEM export extends slightly beyond crop area

tjhb
tjhb Global Mapper UserTrusted User
edited September 2010 in Technical Support
Mike,

I am having trouble cropping a DEM to specified bounds (GM 11.02, build 15 August 2010).

I get an extra half pixel around my area of interest (a total of one extra pixel in each dimension).

The source data is SRTM v4 (from CGIAR) in GeoTIFF tiles, and GTOPO30 (from USGS) in DEM tiles. SRTM tiles are not aligned to round latitude/longitude values.--But my area of interest is, and I need pixel-perfect alignment for my workflow.

Resampling method is set to bicubic for all layers. Output resolution is 2.197265625 seconds per pixel (0.0006103515625 degrees per pixel), which requires resampling up from 3 seconds per pixel (for SRTM) or from 30 seconds (for GTOPO30).

Whether I select an area AOI using the Digitizer tool, and specify "crop to selected area features", or specify export bounds manually (NSWE), the resulting export is not cropped exactly to the specified area, but includes an extra half pixel in each direction. Export boundaries are not aligned with the edges of my AOI, and there is one pixel more in X and in Y than I had calculated.

I would understand this if the option "Export: Keep pixel if any part is inside crop area" were enabled, but it is disabled.

I'm currently trying the same export by script.

Does this ring any bells with you? Is it expected behaviour? Should I try another build, or the current 12 beta?

How can I force GM to resample source DEMs only within the area of interest I specify--so that pixel edges will coincide with the edges of that area?

Regards,

Tim

Comments

  • tjhb
    tjhb Global Mapper User Trusted User
    edited September 2010
    The script export has completed. I get exactly the same result as from the GUI methods.
  • tjhb
    tjhb Global Mapper User Trusted User
    edited September 2010
    I also get the same result (in the GUI) using nearest neighbour resampling on all source layers instead of bicubic.
  • global_mapper
    global_mapper Administrator
    edited September 2010
    Keep in mind that elevation exports take as their bounding box the coordinates for the center of each grid cell, so if your specified export bounds are right on degree boundaries the actual grid cell values will be on those degree boundaries, but when visualizing the data it will expand half a pixel in each direction since the values are at grid cell centers. If you need the outside of the grid cells on even boundaries and not the actual sample locations, you will need to shrink your export bounds by one half sample spacing in each direction.

    Also by default your sample spacing is preserved over your bounds, so your bounds may get expanded by a fraction of a pixel on the right and bottom sides. To have the export bounds be maintained at the expense of the sample spacing, check the option in the Advanced section of the General tab of the Configuration dialog to maintain export bounds over sample spacing.

    Let me know if I can be of further assistance.

    Thanks,

    Mike
    Global Mapper Support
    support@globalmapper.com
  • tjhb
    tjhb Global Mapper User Trusted User
    edited September 2010
    Thanks Mike. That gives me enough to work with.
    Keep in mind that elevation exports take as their bounding box
    the coordinates for the center of each grid cell

    Is that the case regardless of the data model, i.e. even for elevation data where "pixel is area" (as indeed is reported in this case) as well as when "pixel is point"?

    Is there a good reason? It makes it somewhat difficult to accurately combine data supplied in elevation format with data that is nominally imagery. (Quite often the distinction is arbitrary. Landcover data for example can be delivered in either form.)
    If you need the outside of the grid cells on even boundaries and not the actual sample locations,
    you will need to shrink your export bounds by one half sample spacing in each direction.

    I will try this.

    Just upgraded to version 12 by the way, so I will enjoy using that today.
  • global_mapper
    global_mapper Administrator
    edited September 2010
    All gridded layers (whether they actually contain elevation data or just look that way) will treat the data as "pixel is point". If the data set specified "pixel is area" the coordinate space is converted to "pixel is point" on import. So you know if you have a gridded layer that is being treated as terrain, then it's being treated as "pixel is point". Typically for data sets where it is not clear if the data should be treated as imagery or elevation data (i.e. landcover data in say a flat binary file), Global Mapper will ask the user whether the data is an image or terrain.

    Let me know if I can be of further assistance.

    Thanks,

    Mike
    Global Mapper Support
    support@globalmapper.com
  • tjhb
    tjhb Global Mapper User Trusted User
    edited September 2010
    By the way, one thing that would be really handy is to be able to switch, for both elevation data and imagery, between a "pixel is area" interpretation and "pixel is point". And perhaps more important, always to be able to read which interpretation is being used, in the layer's metadata. I'm not certain but think it's currently possible only for those formats that store the information, like GeoTIFF.

    Switching models wouldn't often be required, but sometimes it is, if data is delivered in the wrong format, or to align data that uses different formats. One current way of doing it is to re-register to a position offset one half pixel.
  • tjhb
    tjhb Global Mapper User Trusted User
    edited September 2010
    All gridded layers (whether they actually contain elevation data or just look that way) will treat the data as "pixel is point".
    If the data set specified "pixel is area" the coordinate space is converted to "pixel is point" on import. So you know if you have a gridded layer
    that is being treated as terrain, then it's being treated as "pixel is point".

    Thanks I will try to digest this.

    Would it make sense to have an option to treat gridded data as "pixel is area" (like an image), or would that be too disruptive?

    You have probably thought about this and have a good reason why "pixel is point" is always imposed.

    Why had I never realised any of this? It has implications for working with other software (in my case, mainly Manifold) that assumes the "pixel is area" model for gridded data.
  • tjhb
    tjhb Global Mapper User Trusted User
    edited September 2010
    If the data set specified "pixel is area" the coordinate space is converted to "pixel is point" on import.

    What does that conversion amount to Mike? (I'm full of questions.)
  • global_mapper
    global_mapper Administrator
    edited September 2010
    While it would be technically possible to have an option to treat gridded data sets as "pixel is area", it would have a lot of implications in a lot of places in the code that could make it problematic. It seems odd that Manifold uses "pixel is area" for gridded data sets, I thought they were almost universally treated as "pixel is point" as it seems to make more logical sense to provide the coordinates of the actual grid sample location vs. the edge of the sample cell.

    The conversion from pixel is point to pixel is area is very simple, basically just a half sample spacing shift.

    Thanks,

    Mike
    Global Mapper Support
    support@globalmapper.com
  • tjhb
    tjhb Global Mapper User Trusted User
    edited September 2010
    While it would be technically possible ... it would have a lot of implications
    in a lot of places in the code that could make it problematic.

    Understood, and that's not surprising. Would it be as disruptive to provide an export option to write a DEM as "pixel is area"? That would be very welcome.
    I thought [gridded data sets] were almost universally treated as "pixel is point" as it seems to make more logical sense
    to provide the coordinates of the actual grid sample location vs. the edge of the sample cell.

    It makes sense if the grid was created by sampling that way. But not if it was created by scanning or averaging over a pixel-sized window, as many grids are. The fact that the CGIAR version of SRTM uses "pixel is area" is an obvious, and authoritative, counter-example. In general, where metadata explicitly declares a grid to be "pixel is area", it seems odd to assume (in effect) that this is a mistake.
    The conversion from pixel is point to pixel is area is very simple, basically just a half sample spacing shift.

    Simple, but it's mildly alarming, to learn that Global Mapper is silently repositioning data on import, as if it knows better than the user or the metadata how the data should be interpreted. And it's unhelpful (misleading) to see, for a GeoTIFF DEM that was imported with "pixel is area" registration, metadata still describing it is "pixel is area" if that is no longer the case. This data shifting has important implications for using GM to export gridded data for supply to clients.

    I think these are subjects deserving clear coverage (even a gentle warning) in the user manual.

    I hope there might be some room for change in the behaviour, at least on export, without fundamental disruption to your code base.
  • global_mapper
    global_mapper Administrator
    edited September 2010
    I think there may be some misunderstanding as to how Global Mapper treats "pixel is area" DEMs on import. It doesn't shift the data, it just notes that the coordinates are "pixel is area" and then determines the "pixel is point" coordinate that would correspond to the same placement. Global Mapper is fully aware that different data sets use different "pixel is" coordinate assignments, it just normalizes them all to "pixel is point" on import to get a consistent treatment inside of Global Mapper. For the user it really shouldn't matter how the data is actually stored in the file, Global Mapper automatically handles it and makes sure the data is put in the right place. There should be no actual shift in the placement of the data.

    For GeoTIFF the metadata listed indicates what is in the GeoTIFF file itself. The raw coordinates from the GeoTIFF file will be listed as well (if present), and the "pixel is area" or "pixel is point" indication would apply to those, not to the bounding coordinates listed for the layer which are the same across all gridded types and not specific to GeoTIFF.

    It could be possible to create a "pixel is area" export for some DEM formats (in fact some DEM formats already specifically define the coordinates as the cell corner rather than the cell center, so that is already done for those). Most formats define a specific convention, GeoTIFF is the only one that comes to mind where either method can be used, although it seems that it wouldn't matter much as any compliant GeoTIFF reader would be able to interpret the GeoTIFF header to properly place the file regardless of which convention was used.

    Let me know if I can be of further assistance.

    Thanks,

    Mike
    Global Mapper Support
    support@globalmapper.com
  • tjhb
    tjhb Global Mapper User Trusted User
    edited September 2010
    Thanks for clearing up my confusion—and anxiety. No data shifting: only a switch of reference system. That makes perfect sense. It also makes sense that you'd only have one internal reference system.

    Your second paragraph makes sense too. Thanks for your patient explanations.

    So for present purposes all I'm wanting is to find an export format that uses the "Pixel is Area" standard. (Not GeoTIFF in this case, because the data is too big and Manifold can't read BigTIFF.)

    I had been using ESRI FLT. (Does that use "Pixel is Point" by definition or is there an "Area" option?)

    Do you know of a reference listing DEM formats using "Pixel is Area"? I suspect MapInfo .tab format does; whether it can cope with a large DEM I don't know (I'll try). Any format supported by Manifold would do me. (Those are listed here, under the heading Surface Formats: http://www.manifold.net/info/formats.shtml.)
  • tjhb
    tjhb Global Mapper User Trusted User
    edited September 2010
    Not TAB (since GM doesn't export to this format)...

    Not BIL (a test shows the result is Pixel is Point)...

    Bingo. Looks like Idrisi RST format does the trick, subject to size.
  • tjhb
    tjhb Global Mapper User Trusted User
    edited September 2010
    No problem with file size, and no problem reading the Idrisi file in Manifold.

    But there are problems.

    I specify an export to Idrisi format at 0.0006103515625 degrees per pixel in X and Y, with "Always generate square pixels" on, and the General options "Maintain export bounds instead of sample spacing" and "Export: Keep pixel if any part is within crop area" both off. For my AOI I used a rectangular area with its northwest corner at exactly 3.75°E, 57.5°N and southest corner at exactly 20°E, 33.75°N.

    The target dimensions are therefore 16.25 x 23.75 degrees in X and Y. Dividing by gives 26624 by 38912 pixels (we should probably call them cells, but I'll stick to pixels for consistency).

    (1) Global Mapper exports one too many pixels in X and Y, 26625 and 38913—the same number of pixels it would export for a "Pixel is Point" registration, for example to FLT (with edge cells lying half outside the AOI).

    (2) After importing the GM result, Manifold correctly (perhaps unthinkingly) interprets registration as "Pixel is Area". And the outer edges of all edge pixels are correctly aligned to the AOI used for cropping. But because of point 1, this must be at the expense of pixel resolution, which is reduced proportionately to 0.00061032863849765258 in X and 0.00061033587747025416 in Y.

    (3) Global Mapper does not interpret the result on a "Pixel is Area" basis. (If this is the only option for Idrisi, then it probably should.) After loading the export back in, cell edges are not aligned to the AOI used for cropping, but with a half-cell margin (so apparently as Pixel is Point). Pixel resolution is correctly read, however (or very nearly, within what looks like a binary rounding error: 0.000610351562515262 in X and 0.00061035156251526 in Y).

    To me point 1 looks like a bug. Point 3 may or may not be.

    Maybe this can be fixed. In the meantime, instead of trying to force a "Pixel is Area" export, I'm going back to an earlier suggestion of Mike's:
    If you need the outside of the grid cells on even boundaries and not the actual sample locations,
    you will need to shrink your export bounds by one half sample spacing in each direction.

    Only instead of shrinking the AOI, I'm enlarging it by a buffer of half the pixel spacing (0.00030517578125 degrees).

    This gives me an export (on the Pixel is Point model, to FLT) having an extra two pixels in X and Y (i.e. a margin of one whole pixel). (This is useful because I'll be calculating slope and aspect from the result. The one-pixel margin can easily be trimmed from the resulting images in Photoshop.)

    This works! Net of the one-pixel margin, pixel boundaries are aligned to my original (unbuffered) AOI. Thank heavens for that.

    Sometimes I think life should be easier, though.

    // It strikes me that if Mike wanted to add an option to export any elevation data to a Pixel is Area registration model, it might be enough in practice for him to simply shrink the specified area of interest inwards by half a pixel, in line with the suggested workaround. Or almost enough. I don't know what should happen to data that was not being cropped.
  • tjhb
    tjhb Global Mapper User Trusted User
    edited September 2010
    Above, in paragraph 4, "Dividing by gives..." should read "Dividing by 0.0006103515625 gives...".
  • global_mapper
    global_mapper Administrator
    edited September 2010
    What makes you think that the Idrisi RST format uses "pixel is area"? In Global Mapper both imagery and elevation Idrisi RST files interpret the coordinates as "pixel is point". I can't find any documentation anywhere of how they are supposed to be interpreted, but I must have had a reason to interpret the imagery specifically as "pixel is point" from the coordinates in the RDC file that comes with a RST.

    For any gridded elevation export from Global Mapper the coordinates that you specify are "pixel is point" coordinates, regardless of whether the destination file coordinates will store the values as "pixel is point" or "pixel is area". Otherwise it could get very confusing as users would need to know what coordinate convention the destination file format used, and most users won't know what the different conventions even are (nor should they need to know, software should handle this).

    It would certainly be possible to allow specification of elevation exports as "pixel is area" coordinates by subtracting a half sample spacing for the export, it would just be a matter of where to put the option for "pixel is area" export. I wouldn't want to clutter up the normal export interface as there are probably only a very tiny number of users that would ever use it. Maybe in the Advanced section of the General tab of the Configuration dialog as just a global switch been coordinate references? Does Manifold only support "pixel is area" and not handle either convention?

    Thanks,

    Mike
    Global Mapper Support
    support@globalmapper.com
  • tjhb
    tjhb Global Mapper User Trusted User
    edited September 2010
    What makes you think that the Idrisi RST format uses "pixel is area"?

    Sorry, I seem to have misunderstood again Mike. Your question might be largely rhetorical but I'll answer it anyway. Probably only for two reasons, neither very good, it seems. (1) Because I wanted to find a format for which Pixel is Area was the only option. I drew an inference that some such formats existed from your words
    It could be possible to create a "pixel is area" export for some DEM formats
    (in fact some DEM formats already specifically define the coordinates as the cell corner
    rather than the cell center, so that is already done for those)

    But no doubt my interpretation there was wishful thinking.

    (2) Because of how resulting RST surfaces appear after import to Manifold, compared to the appearance of FLT surfaces. If I import into Manifold an FLT exported from GM, the outer edges of the edge pixels fall outside the AOI used for cropping, by half the sample spacing. If I do the same with an RST exported by GM, then in Manifold the outer edges are aligned exactly to the AOI. (At first I didn't notice that N+1 and M+1 pixels were being scaled into this extent, reducing sample spacing.)

    Bearing in mind your reply, I'm very unsure whether Global Mapper, or Manifold, or both or neither, is interpreting the RST format correctly. They do seem to interpret it differently.

    Since I have almost no knowledge of Idrisi products I can't venture a view, except to say (what I've been thinking all along) that it's no wonder it's so common to read of datasets like SRTM, ETOPO and what not being corrected in successive versions for half-pixel data shifts. This all makes my brain hurt. I'm really grateful for your continued patience in discussing these matters (especially since I think it's getting late where you are). I hope if we get this ironed out and simplified, at least in words, it might help others who find their data doing unexpected things.

    If you added a new option, so much the better (for those who need it). About where to put it:
    Maybe in the Advanced section of the General tab of the Configuration dialog as just a global switch been coordinate references?

    That sounds perfect to me.
    Does Manifold only support "pixel is area" and not handle either convention?

    Yes, only "Pixel is Area" so far, as far as version 8. I'd be surprised if the next version (not far away) didn't allow users to choose which model they want: in effect a choice whether to treat raster data as a grid of points or as a data image (and to change their minds for different purposes). In version 8 this can be done manually, by addressing the data via spatial SQL, where you can define pixel offsets however you like. As you'd expect that tends to be slower than using built-in functions, all of which (including the GPU-enabled Surface Tools) currently reference pixels on the area model. The lower left corner of each pixel is its reference point. A bounding box of a set of cells/pixels includes their combined areal coverage.
  • global_mapper
    global_mapper Administrator
    edited September 2010
    I'm pretty sure there are a few elevation format that use "pixel is area" for referencing, but I can't seem to think of them off the top of my head. I think they were some odd formats so it's quite likely something that Manifold doesn't support.

    I think what I will do is add the option to specify elevation export bounds as "pixel is area", then if that option is checked I will update the elevation GeoTIFF export to use "pixel is area" method since that format I know supports both. You would still be able to specify the export coordinates as "pixel is area" for any elevation export with the option checked, the conversion to the appropriate reference method would just be done automatically for you when the export is done based on what the format being export to specifies should be used. I'll let you know when I have this ready for test.

    Thanks,

    Mike
    Global Mapper Support
    support@globalmapper.com
  • global_mapper
    global_mapper Administrator
    edited September 2010
    I have added a new option in the Advanced section of the General tab of the Configuration dialog (about the 5th option down) to indicate that elevation export bounds are specified as "pixel is area" rather than "pixel is point". This is available in a new v12 beta installer available for testing at:

    32-bit: http://www.globalmapper.com/global_mapper12_setup.exe
    64-bit: http://www.globalmapper.com/global_mapper12_setup_64bit.exe

    Let me know if I can be of further assistance.

    Thanks,

    Mike
    Global Mapper Support
    support@globalmapper.com
  • tjhb
    tjhb Global Mapper User Trusted User
    edited September 2010
    What a nice surprise to wake up to. How did you do that so quickly? I'll test it over the next few days with all the formats we normally use and let you know of any problems.

    This option makes a real difference to us for interoperability, not only with other software such as Manifold and Photoshop but also for different client purposes. Now I know to ask myself (but only occasionally a client, since most won't know) what method of pixel registration is appropriate to a given data set and purpose, and we can adjust very easily. Previously I've just assumed the problem away.

    Thanks Mike (loudly).