Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I have some data looking like this:

<?xml version="1.0" encoding="utf-8"?>
<levels>
  <level nr="0" tileWidth="512" tileHeight="512" envSizeX="5000" envSizeY="5000" minX="465000" maxX="470000" minY="105000" maxY="110000" scaleFactor="1" tilesCountX="10" tilesCountY="10" />
  <level nr="1" tileWidth="512" tileHeight="512" envSizeX="4613" envSizeY="4613" minX="465000" maxX="469999" minY="105001" maxY="110000" scaleFactor="1.0837333333333332" tilesCountX="10" tilesCountY="10" />
  <level nr="2" tileWidth="512" tileHeight="512" envSizeX="2286" envSizeY="2286" minX="465000" maxX="469955" minY="105045" maxY="110000" scaleFactor="2.1674666666666664" tilesCountX="5" tilesCountY="5" />
  <level nr="3" tileWidth="512" tileHeight="512" envSizeX="1134" envSizeY="1134" minX="465000" maxX="469915" minY="105085" maxY="110000" scaleFactor="4.3349333333333329" tilesCountX="3" tilesCountY="3" />
  <level nr="4" tileWidth="512" tileHeight="512" envSizeX="563" envSizeY="563" minX="465000" maxX="469879" minY="105121" maxY="110000" scaleFactor="8.6698666666666657" tilesCountX="2" tilesCountY="2" />
  <level nr="5" tileWidth="512" tileHeight="512" envSizeX="275" envSizeY="275" minX="465000" maxX="469777" minY="105223" maxY="110000" scaleFactor="17.339733333333331" tilesCountX="1" tilesCountY="1" />
  <level nr="6" tileWidth="512" tileHeight="512" envSizeX="134" envSizeY="134" minX="465000" maxX="469658" minY="105342" maxY="110000" scaleFactor="34.679466666666663" tilesCountX="1" tilesCountY="1" />
  <level nr="7" tileWidth="512" tileHeight="512" envSizeX="67" envSizeY="67" minX="465000" maxX="469623" minY="105377" maxY="110000" scaleFactor="69.358933333333326" tilesCountX="1" tilesCountY="1" />
  <level nr="8" tileWidth="512" tileHeight="512" envSizeX="33" envSizeY="33" minX="465000" maxX="469554" minY="105446" maxY="110000" scaleFactor="138.71786666666665" tilesCountX="1" tilesCountY="1" />
  <level nr="9" tileWidth="512" tileHeight="512" envSizeX="16" envSizeY="16" minX="465000" maxX="469432" minY="105568" maxY="110000" scaleFactor="277.4357333333333" tilesCountX="1" tilesCountY="1" />
  <level nr="10" tileWidth="512" tileHeight="512" envSizeX="8" envSizeY="8" minX="465000" maxX="469432" minY="105568" maxY="110000" scaleFactor="554.87146666666661" tilesCountX="1" tilesCountY="1" />
</levels>

I have to use this data to designate tiles' extents (minX, maxX, minY, maxY). I wrote a method that gets a level and calculates its tiles. Then it gets the tiles and chooses these which belong to a viewport area.

I thought that everything works fine, but I noticed the differences between data (when I compare GdalInfo data from metadata of TIFF file and my counted values there are some differences...

e.g.
GDAL_INFO = 465000, 465554.871466667, 109445.128533333, 110000
MY_EXTENTS = (465000, 465499.9, 105001, 105500.9)

Here's my method:

public static IEnumerable<Envelope> GetTiles(string xmlFileName, int level, Rect viewportData, bool debug = false)
{
    // viewport data
    double x = viewportData.X, y = viewportData.Y;
    double width = viewportData.Width, height = viewportData.Height;

    // load level data
    XDocument doc = XDocument.Load(xmlFileName);
    if (doc == null)
    {
        throw new FileNotFoundException("File '" + xmlFileName + "' not found.");
    }

    XElement levelData =
        doc.Element("levels")
            .Descendants()
            .First(element => int.Parse(element.Attribute("nr").Value) == level);

    double minX = double.Parse(levelData.Attribute("minX").Value);
    double maxX = double.Parse(levelData.Attribute("maxX").Value);
    double minY = double.Parse(levelData.Attribute("minY").Value);
    double maxY = double.Parse(levelData.Attribute("maxY").Value);

    int tilesCountX = int.Parse(levelData.Attribute("tilesCountX").Value);
    int tilesCountY = int.Parse(levelData.Attribute("tilesCountY").Value);

    Envelope[][] extents = new Envelope[tilesCountX][];

    Envelope bounds = new Envelope();
    bounds.MinX = 0;
    bounds.MaxX = 0;
    bounds.MinY = 0;
    bounds.MaxY = 0;
    for (int i = 0; i < extents.Length; i++)
    {
        extents[i] = new Envelope[tilesCountY];
        for (int j = 0; j < extents[i].Length; j++)
        {
            // count extents for pieces
            double tileGeoWidth = (maxX - minX)/(double) tilesCountX;
            double tileGeoHeight = (maxY - minY)/(double) tilesCountY;
            double offsetX = (double) j*tileGeoWidth;
            double offsetY = (double) i*tileGeoHeight;

            extents[i][j] = new Envelope();
            extents[i][j].MinX = minX + offsetX;
            extents[i][j].MaxX = minX + tileGeoWidth + offsetX; // poprawic dla ostatnich (niewymiarowych) tilesow
            extents[i][j].MinY = minY + offsetY;
            extents[i][j].MaxY = minY + tileGeoHeight + offsetY;

            // check if a piece matches to the viewport
            // 1. find geo-bounds
            if (offsetX <= x)
                bounds.MinX = offsetX;
            if (offsetX <= x + width)
                bounds.MaxX = offsetX;
            if (offsetY <= y)
                bounds.MinY = offsetY;
            if (offsetY <= y + height)
                bounds.MaxY = offsetY;
        }
    }

    // 2. increase "max" bounds (+ width/height)
    bounds.MaxX += width;
    bounds.MaxY += height;

    var intersectedTiles =
        extents
        .SelectMany(es => es)
        .Where(e => EnvIntersects(e, bounds))
        .ToList();

    if (debug)
    {
        Console.WriteLine("GLOBAL_EXTENT = {0}, {1}, {2}, {3}", minX, maxX, minY, maxY);
        Console.WriteLine("TILES_COUNT = {0}/{1}", tilesCountX, tilesCountY);
        for (int i = 0; i < extents.LongLength; i++)
        {
            for (int j = 0; j < extents.Length; j++)
            {
                Console.WriteLine("EXTENTS = ({0}, {1}, {2}, {3})", extents[i][j].MinX, extents[i][j].MaxX, extents[i][j].MinY, extents[i][j].MaxY);
            }
        }
        foreach (Envelope intersectedTile in intersectedTiles)
        {
            Console.WriteLine("INTERSECTION = ({0}, {1}, {2}, {3})", intersectedTile.MinX, intersectedTile.MaxX,
                              intersectedTile.MinY, intersectedTile.MaxY);
        }
    }

    return intersectedTiles;
}

And an example call:

Rect viewportData = new Rect(100000, 200000, 70000, 50000); 
GdalRetile.GetTiles(@"D:\#PYRAMID\su60ne.xml", 1, viewportData, true);

Could anyone help me to correct this code?
Thanks in advance!

share|improve this question

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Browse other questions tagged or ask your own question.