[gdal-dev] Re: How to divide a raster layer into a set of tiles

[Eli Adam]

Gilles,
     I do a similar process (build large vrt mosaic, reproject the
mosaic vrt to another vrt, tile the reprojected mosaic vrt in new
projection to tif) although my end tiles are in tif and not compressed. 
My method is the brute force way and not elegant at all.  You may also
want to know that I have been advised
(http://news.gmane.org/gmane.comp.gis.gdal.devel) that this is very
inefficient  and pushing the VRT format.  Steps 1 and 2 are the same as
you describe.

1) create a single GDAL layer from all initial tiles with gdalbuildvrt
2) reproject this virtual raster with gdalwarp (outputs a new .vrt)
3) Tile it with gdal_translate:  gdal_translate -projwin 7288440 446648
7306569 429923 mosaicvrtprj.vrt tile1.tif

The brute force part is between 2 and 3.  Make a vector file (shp or
PostGIS) grid of the new tiles in the desired projection.  Convert the
polygon grid to points, then using SQL in a database, group by the
polygon name (which needs to be retained when converted to points) and
take the min/max of x and y.  Then use SQL to concatenate gdal_translate
commands with the correct parameters.  The end result is a table of
gdal_translate commands that will make the tiles.  Save as text and make
needed minimal modification to fire off all those commands.   Not an
elegant method at all.  It could probably be done all in bash with
gdalinfo, grep, and sed, which would at least not require a database,
even if it still is a brute force method.  

I don't know if using ecw and other formats or compression options may
prevent the reading through of vrt files.  My work has been simpler,
uncompressed tif from AI binary grids.  You may also want to make sure
that each of your vrt files has all the proper projection and other
information from steps 1 and 2.  gdalinfo ortho-dept04-2004-lam93.vrt
might help.  


Bests, Eli


   2. Re: How to divide a raster layer into a set of tiles?
      (Gilles Bassi?re)
   3. Re: Re: How to divide a raster layer into a set of tiles?
      (Gilles Bassi?re)
   4. Re: How to divide a raster layer into a set of tiles?
      (Jukka Rahkonen)

------------------------------

Message: 2
Date: Tue, 23 Mar 2010 16:01:33 +0100
From: Gilles Bassi?re <gbassiere at gmail.com>
Subject: [gdal-dev] Re: How to divide a raster layer into a set of
	tiles?
To: Christian M?ller <christian.mueller at nvoe.at>
Cc: gdal-dev at lists.osgeo.org 
Message-ID: <4BA8D7CD.1050300 at gmail.com>
Content-Type: text/plain; charset=ISO-8859-1

Christian Müller wrote:
> That is a surprise.
> In all my tests I had no memory problems with gdal_retile. Sometimes
I
> had some memory problems with gdal_merge. The biggest input file for
> gdal_retile was a 80 GB Erdas image (not compressed) and it works
fine.
> I have never testet with an ECW File and my tests run on a Linux Box
or
> AIX Box.
> Perhaps the problem is in the ECW gdal driver. You can try to
transform
> your ecw file in another format using gdal_transform.
> 

That's may be caused by something wrong with the ECW driver indeed. I
will test gdal_retile again once I have a clean GDAL+ECW install.

Regards
Gilles

> 
> 
> Gilles Bassière writes:
>> Christian Müller wrote:
>>> gdal_retile.py definitely creates georeferenced tiles, I developed
this
>>> utility.
>>> I did some experiments reprojecting the single tiles while
retiling, but
>>> it simply does not work. As a consequence I canceled the
reprojecting
>>> support.
>>> The solution for your problem is:
>>> 1) Use gdal_merge to create the big image
>>> 2) reproject with gdal_warp
>>> 3) use gdal_retile to create your new tiles
>>> I did a lot of experiments to avoid creating the big picture
without
>>> success.
>>
>> I understand. In my case, gdal_retile.py eventually crashed with a
>> Python MemoryError. I assume it is unsurprising with a 7.5GB input
for a
>> single Python process.
>>>
>>> Even Rouault writes:
>>>> Le Monday 22 March 2010 19:35:47 Gilles Bassière, vous avez écrit
:
>>>>> Hi,
>>>>> The problem with this workflow is that it creates
non-georeferenced
>>>>> tiles (I don't know if this is the intended behaviour of this
>>>>> command).
>>>>
>>>> gdal_retile.py *does* create georeferenced tiles. I've just
verified
>>>> it. What make you think the contrary ?
>>>> _______________________________________________
>>>> gdal-dev mailing list
>>>> gdal-dev at lists.osgeo.org 
>>>> http://lists.osgeo.org/mailman/listinfo/gdal-dev 
>>>  
>>>
>>> _______________________________________________
>>> gdal-dev mailing list
>>> gdal-dev at lists.osgeo.org 
>>> http://lists.osgeo.org/mailman/listinfo/gdal-dev 
>>  
>>
>> -- 
>> Gilles Bassière - Web/GIS software engineer
>> http://gbassiere.free.fr/ 
> 
> 


-- 
Gilles Bassière - Web/GIS software engineer
http://gbassiere.free.fr/ 


------------------------------

Message: 3
Date: Tue, 23 Mar 2010 16:08:58 +0100
From: Gilles Bassi?re <gbassiere at gmail.com>
Subject: Re: [gdal-dev] Re: How to divide a raster layer into a set of
	tiles?
To: Jukka Rahkonen <jukka.rahkonen at mmmtike.fi>
Cc: gdal-dev at lists.osgeo.org 
Message-ID: <4BA8D98A.2020905 at gmail.com>
Content-Type: text/plain; charset=UTF-8

Jukka Rahkonen wrote:
> Gilles BassiÃ̈re <gbassiere <at> gmail.com> writes:
> 
>> Hi,
>>
>> I'm trying to reproject an aerial imagery dataset with GDAL.
>>
>> The initial dataset is composed of ECW tiles in EPSG:27572 and my
goal
>> is to compute a new set set of ECW tiles projected in EPSG:2154.
> ...
> 
>> My current workflow consists in the following steps:
>> 1) create a single GDAL layer from all initial tiles with
gdalbuildvrt
>> 2) reproject this virtual raster with gdalwarp (outputs a new .vrt)
> 
> I would perhaps calculate the extents I want the output tiles to have
and run
> gdalwarp for each tile with -te option.
> 
> -Jukka Rahkonen-
> 

Hi Jukka,

Thanks for your suggestion. I've been able to get the extent
automatically computed by a little Bash scripting:

 for x in `seq -w 895 5 1020`; do
     for y in `seq -w 6290 5 6410`; do
         e="$((10#${x}*1000)) $((10#${y}*1000)) $(((10#${x}+5)*1000))
$(((10#${y}+5)*1000))"
         echo "=== Processing tile $x-$y ==="
         gdalwarp -s_srs "+init=epsg:27572" -t_srs "+init=epsg:2154"
-te
$e -ts 10000 10000 ortho-dept04-2004-la2e.vrt
ortho-dept04-2004-lam93-tif/04-2004-$x-$y-lam93.tif
     done
 done

Regards

-- 
Gilles BassiÃ̈re - Web/GIS software engineer
http://gbassiere.free.fr/ 


------------------------------

Message: 4
Date: Tue, 23 Mar 2010 11:54:50 +0000 (UTC)
From: Jukka Rahkonen <jukka.rahkonen at mmmtike.fi>
Subject: [gdal-dev] Re: How to divide a raster layer into a set of
	tiles?
To: gdal-dev at lists.osgeo.org 
Message-ID: <loom.20100323T125053-679 at post.gmane.org>
Content-Type: text/plain; charset=utf-8

Gilles BassiÃ̈re <gbassiere <at> gmail.com> writes:

> 
> Hi,
> 
> I'm trying to reproject an aerial imagery dataset with GDAL.
> 
> The initial dataset is composed of ECW tiles in EPSG:27572 and my
goal
> is to compute a new set set of ECW tiles projected in EPSG:2154.

Hi,

May I ask if you are ready with your new tiles? The task is for sure
doable with
a few alternative ways. It would be nice to know which route you
decided to
follow, perhaps with some examples of the commands you have used.

-Jukka Rahkonen-