[gdal-dev] VRT pixel functions on tiled vrt

classic Classic list List threaded Threaded
9 messages Options
Reply | Threaded
Open this post in threaded view
|

[gdal-dev] VRT pixel functions on tiled vrt

William Kyngesburye
I'm trying to add a pixel function (Python) to a vrt where there are multiple rasters (ComplexSource) tiled together (that is, not stacked like an RGB image).  When I translate the vrt to a tif, only the first ComplexSource is processed, the rest on the raster is set to 0.

I got the basics of the functin from the vrt docs, and it does work, for the first raster tile.  Maybe the python notation for calculating the numpy array is wrong?

The original vrt is created from gdalbuildvrt.  This is the vrt thru the first 2 rasters:

<VRTDataset rasterXSize="14401" rasterYSize="4801">
 <SRS>GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]]</SRS>
 <GeoTransform> -2.5000416666666666e+01,  8.3333333333333263e-04,  0.0000000000000000e+00,  6.7000416666666666e+01,  0.0000000000000000e+00, -8.3333333333333263e-04</GeoTransform>
 <VRTRasterBand dataType="Byte" band="1" subClass="VRTDerivedRasterBand">
   <PixelFunctionType>land</PixelFunctionType>
   <PixelFunctionLanguage>Python</PixelFunctionLanguage>
   <PixelFunctionCode><![CDATA[
import numpy
def land(in_ar, out_ar, xoff, yoff, xsize, ysize, raster_xsize, raster_ysize, buf_radius, gt, **kwargs):
        out_ar[:] = (in_ar[0] > 0)
]]>
   </PixelFunctionCode>
   <ColorInterp>Gray</ColorInterp>
   <ComplexSource>
     <SourceFilename relativeToVRT="1">n63w017.tif</SourceFilename>
     <SourceBand>1</SourceBand>
     <SourceProperties RasterXSize="1201" RasterYSize="1201" DataType="Int16" BlockXSize="1201" BlockYSize="3" />
     <SrcRect xOff="0" yOff="0" xSize="1201" ySize="1201" />
     <DstRect xOff="9600" yOff="3600" xSize="1201" ySize="1201" />
     <NODATA>-32768</NODATA>
   </ComplexSource>
   <ComplexSource>
     <SourceFilename relativeToVRT="1">n63w018.tif</SourceFilename>
     <SourceBand>1</SourceBand>
     <SourceProperties RasterXSize="1201" RasterYSize="1201" DataType="Int16" BlockXSize="1201" BlockYSize="3" />
     <SrcRect xOff="0" yOff="0" xSize="1201" ySize="1201" />
     <DstRect xOff="8400" yOff="3600" xSize="1201" ySize="1201" />
     <NODATA>-32768</NODATA>
   </ComplexSource>
...
 </VRTRasterBand>
</VRTDataset>


-----
William Kyngesburye <kyngchaos*at*kyngchaos*dot*com>
http://www.kyngchaos.com/

The equator is so long, it could encircle the earth completely once.

_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: VRT pixel functions on tiled vrt

Even Rouault-2

On jeudi 26 avril 2018 08:22:05 CEST William Kyngesburye wrote:

> I'm trying to add a pixel function (Python) to a vrt where there are

> multiple rasters (ComplexSource) tiled together (that is, not stacked like

> an RGB image). When I translate the vrt to a tif, only the first

> ComplexSource is processed, the rest on the raster is set to 0.

>

> I got the basics of the functin from the vrt docs, and it does work, for the

> first raster tile. Maybe the python notation for calculating the numpy

> array is wrong?

>

> The original vrt is created from gdalbuildvrt. This is the vrt thru the

> first 2 rasters:

 

William,

 

in a VRTDerivedRasterBand, the sources are considered to be a stack, and in_ar[0] will return the first source and in_ar[1] the seconded one.

 

What you need to do is create a regular tiled VRT, and use that one has the single source of a VRTDerivedRasterBand

 

Even

 

--

Spatialys - Geospatial professional services

http://www.spatialys.com


_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: VRT pixel functions on tiled vrt

William Kyngesburye
So, a vrt of a vrt?

I wanted to avoid that, because as I understand the source, it reads each raster in the derived band as a single numpy array, which would be all the tiled rasters together.  Most of the vrts that I want to process this way are huge, GBs.

I guess I can make a derived band vrt of each tile, then merge those together as a normal vrt.

I already tried running gdal_calc.py on the normal vrts to do the same thing into a geotiff, but it takes forever, while the pixel functions (also python) appear to be much faster.

On Apr 26, 2018, at 12:41 PM, Even Rouault <[hidden email]> wrote:

On jeudi 26 avril 2018 08:22:05 CEST William Kyngesburye wrote:
> I'm trying to add a pixel function (Python) to a vrt where there are
> multiple rasters (ComplexSource) tiled together (that is, not stacked like
> an RGB image). When I translate the vrt to a tif, only the first
> ComplexSource is processed, the rest on the raster is set to 0.
>
> I got the basics of the functin from the vrt docs, and it does work, for the
> first raster tile. Maybe the python notation for calculating the numpy
> array is wrong?
>
> The original vrt is created from gdalbuildvrt. This is the vrt thru the
> first 2 rasters:

 

William,

 

in a VRTDerivedRasterBand, the sources are considered to be a stack, and in_ar[0] will return the first source and in_ar[1] the seconded one.

 

What you need to do is create a regular tiled VRT, and use that one has the single source of a VRTDerivedRasterBand

 

Even

 

--
Spatialys - Geospatial professional services
_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev

-----
William Kyngesburye <kyngchaos*at*kyngchaos*dot*com>

"Time is an illusion - lunchtime doubly so."

- Ford Prefect



_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: VRT pixel functions on tiled vrt

Even Rouault-2

On jeudi 26 avril 2018 12:56:10 CEST William Kyngesburye wrote:

> So, a vrt of a vrt?

 

Yes

 

>

> I wanted to avoid that, because as I understand the source, it reads each

> raster in the derived band as a single numpy array, which would be all the

> tiled rasters together.

 

No, if you use gdal_translate, it will process the source by reasonable chunks, so you should be able to process arbitrarily large sources.

 

Even

 

--

Spatialys - Geospatial professional services

http://www.spatialys.com


_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: VRT pixel functions on tiled vrt

William Kyngesburye
So, if the gdal library queries for a region of the vrt (ie processing it in chunks), the numpy array will be just that region?  I'll give it a try.

What about something like viewing the whole vrt in QGIS - will it be read in chunks as well?  IOW, is it the GDAL library that breaks it into chunks at reading, or is it the responsibility of programs (gdal_translate, QGIS,...) to ask for chunks?

On Apr 26, 2018, at 12:59 PM, Even Rouault <[hidden email]> wrote:

On jeudi 26 avril 2018 12:56:10 CEST William Kyngesburye wrote:
> So, a vrt of a vrt?

 

Yes

 

>
> I wanted to avoid that, because as I understand the source, it reads each
> raster in the derived band as a single numpy array, which would be all the
> tiled rasters together.

 

No, if you use gdal_translate, it will process the source by reasonable chunks, so you should be able to process arbitrarily large sources.

 

Even

 

--
Spatialys - Geospatial professional services
_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev

-----
William Kyngesburye <kyngchaos*at*kyngchaos*dot*com>

"We can die but once, and that once we must die.  To be always fearing, then, would not avert it, and would make life miserable."

- Tarzan, on death


_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: VRT pixel functions on tiled vrt

William Kyngesburye
OK, it worked as vrt of a vrt.  Memory not an issue.

But I have a new problem.  One vrt I want to instead do > -29 (caspian sea, -29 in SRTM data).  The whole region evaluates to 1 (true).  It looks like the in_ar[] is getting truncated to an unsigned byte, 0-255, from its orignal Int16.  Testing some values, anything up to 254 gives me 0/1 areas matching the elevation values in the original DEM, but >255 returns 0 for the whole area (that is, everything >255 was truncated to 255, so nothing is >255).  Thus, > -29 is true for the whole area.

I tried using the python cmp() function, in case the () test was casting each value to Byte, but it complained that it's ambiguous for an array.  The error suggested to use a.any() or a.all(), but I don't know how to do that one.

Here is the script again:

import numpy
def land(in_ar, out_ar, xoff, yoff, xsize, ysize, raster_xsize, raster_ysize, buf_radius, gt, **kwargs):
out_ar[:] = (in_ar[0] > -29)

On Apr 26, 2018, at 1:07 PM, William Kyngesburye <[hidden email]> wrote:

So, if the gdal library queries for a region of the vrt (ie processing it in chunks), the numpy array will be just that region?  I'll give it a try.

What about something like viewing the whole vrt in QGIS - will it be read in chunks as well?  IOW, is it the GDAL library that breaks it into chunks at reading, or is it the responsibility of programs (gdal_translate, QGIS,...) to ask for chunks?

On Apr 26, 2018, at 12:59 PM, Even Rouault <[hidden email]> wrote:

On jeudi 26 avril 2018 12:56:10 CEST William Kyngesburye wrote:
> So, a vrt of a vrt?

 

Yes

 

>
> I wanted to avoid that, because as I understand the source, it reads each
> raster in the derived band as a single numpy array, which would be all the
> tiled rasters together.

 

No, if you use gdal_translate, it will process the source by reasonable chunks, so you should be able to process arbitrarily large sources.

 

Even

 

--
Spatialys - Geospatial professional services
_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev

-----
William Kyngesburye <kyngchaos*at*kyngchaos*dot*com>

"We can die but once, and that once we must die.  To be always fearing, then, would not avert it, and would make life miserable."

- Tarzan, on death

_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev

-----
William Kyngesburye <kyngchaos*at*kyngchaos*dot*com>

First Pogril: Why is life like sticking your head in a bucket filled with hyena offal?
Second Pogril: I don't know.  Why IS life like sticking your head in a bucket filled with hyena offal?
First Pogril: I don't know either.  Wretched, isn't it?

-HitchHiker's Guide to the Galaxy



_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: VRT pixel functions on tiled vrt

William Kyngesburye
I tried:

out_ar[:] = ((in_ar[0] + 29) > 0)

to try to remove the suspected truncation - bring -29 up to 0 before it's truncated - but then elevation 227 is false and everything else is true (<=226, >=228)

> On Apr 27, 2018, at 12:42 PM, William Kyngesburye <[hidden email]> wrote:
>
> OK, it worked as vrt of a vrt.  Memory not an issue.
>
> But I have a new problem.  One vrt I want to instead do > -29 (caspian sea, -29 in SRTM data).  The whole region evaluates to 1 (true).  It looks like the in_ar[] is getting truncated to an unsigned byte, 0-255, from its orignal Int16.  Testing some values, anything up to 254 gives me 0/1 areas matching the elevation values in the original DEM, but >255 returns 0 for the whole area (that is, everything >255 was truncated to 255, so nothing is >255).  Thus, > -29 is true for the whole area.
>
> I tried using the python cmp() function, in case the () test was casting each value to Byte, but it complained that it's ambiguous for an array.  The error suggested to use a.any() or a.all(), but I don't know how to do that one.
>
> Here is the script again:
>
> import numpy
> def land(in_ar, out_ar, xoff, yoff, xsize, ysize, raster_xsize, raster_ysize, buf_radius, gt, **kwargs):
> out_ar[:] = (in_ar[0] > -29)
>
>> On Apr 26, 2018, at 1:07 PM, William Kyngesburye <[hidden email]> wrote:
>>
>> So, if the gdal library queries for a region of the vrt (ie processing it in chunks), the numpy array will be just that region?  I'll give it a try.
>>
>> What about something like viewing the whole vrt in QGIS - will it be read in chunks as well?  IOW, is it the GDAL library that breaks it into chunks at reading, or is it the responsibility of programs (gdal_translate, QGIS,...) to ask for chunks?
>>
>>> On Apr 26, 2018, at 12:59 PM, Even Rouault <[hidden email]> wrote:
>>>
>>> On jeudi 26 avril 2018 12:56:10 CEST William Kyngesburye wrote:
>>> > So, a vrt of a vrt?
>>>  
>>> Yes
>>>  
>>> >
>>> > I wanted to avoid that, because as I understand the source, it reads each
>>> > raster in the derived band as a single numpy array, which would be all the
>>> > tiled rasters together.  
>>>  
>>> No, if you use gdal_translate, it will process the source by reasonable chunks, so you should be able to process arbitrarily large sources.
>>>  
>>> Even
>>>  
>>> --
>>> Spatialys - Geospatial professional services
>>> http://www.spatialys.com
>>> _______________________________________________
>>> gdal-dev mailing list
>>> [hidden email]
>>> https://lists.osgeo.org/mailman/listinfo/gdal-dev
>>
>> -----
>> William Kyngesburye <kyngchaos*at*kyngchaos*dot*com>
>> http://www.kyngchaos.com/
>>
>> "We can die but once, and that once we must die.  To be always fearing, then, would not avert it, and would make life miserable."
>>
>> - Tarzan, on death
>>
>> _______________________________________________
>> gdal-dev mailing list
>> [hidden email]
>> https://lists.osgeo.org/mailman/listinfo/gdal-dev
>
> -----
> William Kyngesburye <kyngchaos*at*kyngchaos*dot*com>
> http://www.kyngchaos.com/
>
> First Pogril: Why is life like sticking your head in a bucket filled with hyena offal?
> Second Pogril: I don't know.  Why IS life like sticking your head in a bucket filled with hyena offal?
> First Pogril: I don't know either.  Wretched, isn't it?
>
> -HitchHiker's Guide to the Galaxy
>
>
> _______________________________________________
> gdal-dev mailing list
> [hidden email]
> https://lists.osgeo.org/mailman/listinfo/gdal-dev

-----
William Kyngesburye <kyngchaos*at*kyngchaos*dot*com>
http://www.kyngchaos.com/

[Trillian]  What are you supposed to do WITH a maniacally depressed robot?

[Marvin]  You think you have problems?  What are you supposed to do if you ARE a maniacally depressed robot?  No, don't try and answer, I'm 50,000 times more intelligent than you and even I don't know the answer...

- HitchHiker's Guide to the Galaxy


_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: VRT pixel functions on tiled vrt

Even Rouault-2
In reply to this post by William Kyngesburye

On vendredi 27 avril 2018 12:42:16 CEST William Kyngesburye wrote:

> OK, it worked as vrt of a vrt. Memory not an issue.

>

> But I have a new problem. One vrt I want to instead do > -29 (caspian sea,

> -29 in SRTM data). The whole region evaluates to 1 (true). It looks like

> the in_ar[] is getting truncated to an unsigned byte, 0-255, from its

> orignal Int16.

 

Check the various data types mentionned in the VRT. Perhaps there's one set to Byte.

It would also probably be safer to define

<SourceTransferType>GDT_Int16</SourceTransferType>

at the same level than the <PixelFunction...> elements

 

--

Spatialys - Geospatial professional services

http://www.spatialys.com


_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: VRT pixel functions on tiled vrt

William Kyngesburye
I missed the SourceTransferType parameter, that worked if I used a simple Int16.

thanks

> On Apr 27, 2018, at 2:32 PM, Even Rouault <[hidden email]> wrote:
>
> On vendredi 27 avril 2018 12:42:16 CEST William Kyngesburye wrote:
> > OK, it worked as vrt of a vrt.  Memory not an issue.
> >
> > But I have a new problem.  One vrt I want to instead do > -29 (caspian sea,
> > -29 in SRTM data).  The whole region evaluates to 1 (true).  It looks like
> > the in_ar[] is getting truncated to an unsigned byte, 0-255, from its
> > orignal Int16.  
>  
> Check the various data types mentionned in the VRT. Perhaps there's one set to Byte.
> It would also probably be safer to define
> <SourceTransferType>GDT_Int16</SourceTransferType>
> at the same level than the <PixelFunction...> elements
>  
> --
> Spatialys - Geospatial professional services
> http://www.spatialys.com
> _______________________________________________
> gdal-dev mailing list
> [hidden email]
> https://lists.osgeo.org/mailman/listinfo/gdal-dev

-----
William Kyngesburye <kyngchaos*at*kyngchaos*dot*com>
http://www.kyngchaos.com/

First Pogril: Why is life like sticking your head in a bucket filled with hyena offal?
Second Pogril: I don't know.  Why IS life like sticking your head in a bucket filled with hyena offal?
First Pogril: I don't know either.  Wretched, isn't it?

-HitchHiker's Guide to the Galaxy


_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev