[gdal-dev] Re: How to divide a raster layer into a set of tiles
Christian Müller
christian.mueller at nvoe.at
Wed Mar 24 01:40:40 EDT 2010
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