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

Eli Adam EAdam at co.lincoln.or.us
Wed Mar 24 10:36:55 EDT 2010


No good reason.  I'm more familiar with the other commands since I use
them often is the main reason. A side benefit is that I can match
existing tiles on other datasets and adopt existing tile names too.  Now
that I've read the utility pages for gdal_merge and gdal_retile I see
that would be an easier method.  I'll try that next time and report if I
have problems.  The reason I recommended it to Gilles was that the
method he was trying was not going smoothly. 

Thanks for contributing these utilities.  

Bests, Eli

>>> "Christian Müller" <christian.mueller at nvoe.at> 3/23/2010 10:40 PM
>>>
A short question 

Why did you not use 

1) gdal_merge
2) gdal_translate
3) gdal_retile 

Btw, if there is no need for reprojection, gdal_retile alone does the
job. 

Did you have problems using these utilities ? 

Eli Adam writes: 

> 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- 
> 
>  
> 
> 
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org 
> http://lists.osgeo.org/mailman/listinfo/gdal-dev 
 



More information about the gdal-dev mailing list