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
add comment

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.