NDVI Calculation from two bands within one Raster

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

NDVI Calculation from two bands within one Raster

JamesH
I am attempting to complete a NDVI calculation using two bands (3 and 4) from a Multiband Raster stored as a Table (nclraster1).

From a previous question asked I am basing my calculation on:

SELECT
        ST_MapAlgebra(ST_MapAlgebra(a.rast, 4, b.rast, 3, "rast - b.rast", "8BUI"), ST_MapAlgebra(a.rast, 4, b.rast, 3, "a.rast + b.rast2", "8BUI"), "a.rast / b.rast", "8BUI")
FROM
        nclraster1 a, nclraster1 b;

This does not work for me, returning the column a.rast does not exist.

Any pointers on where I'm going wrong much appreciated.

James

GIS Undergraduate
Reply | Threaded
Open this post in threaded view
|

Re: NDVI Calculation from two bands within one Raster

Pierre Racine-2
> SELECT
> ST_MapAlgebra(ST_MapAlgebra(a.rast, 4, b.rast, 3, "rast - b.rast",
> "8BUI"), ST_MapAlgebra(a.rast, 4, b.rast, 3, "a.rast + b.rast2", "8BUI"), "a.rast /
> b.rast", "8BUI") FROM
> nclraster1 a, nclraster1 b;
>
> This does not work for me, returning the column a.rast does not exist.

1) Do you have a column of type "raster" in table "nclraster1"?

2) Is it named "rast"?

3) Is this raster tiled? (or Do you have more than one row in "nclraster1"?)

4) What if the results of

SELECT (ST_SummaryStats(rast)).*
FROM nclraster1

?

5) What is the result of:

SELECT PostGIS_Full_Version();

6) Your query should probably more look like (if you have a recent version of PostGIS):

SELECT ST_MapAlgebraExpr(a.rast, 4, b.rast, 3, '([rast1] - [rast2]) / ([rast1] + [rast2])::float')
FROM nclraster1 a, nclraster1 b;

Pierre
_______________________________________________
postgis-users mailing list
[hidden email]
http://postgis.refractions.net/mailman/listinfo/postgis-users
Reply | Threaded
Open this post in threaded view
|

Re: NDVI Calculation from two bands within one Raster

JamesH
> 1) 2) Do you have a column of type "raster" in table "nclraster1"?
Yes rid and rast.

> 3) Is this raster tiled? (or Do you have more than one row in "nclraster1"?)
No , only 1 row in nclraster1

> 4) What if the results of SELECT (ST_SummaryStats(rast)).*FROM nclraster1
count(bigint)= 49729, sum(dbprecision)= 432970, mean(dbprecision)= 8.706589..., stdev(dbprecision)= 10.15282..., min(dbprecision)= 0, max(dbprecision)= 236

> 5) What is the result of:SELECT PostGIS_Full_Version();
"POSTGIS="2.0.0SVN" GEOS="3.3.2dev-CAPI-1.7.2" PROJ="Rel. 4.6.1, 21 August 2008" GDAL="GDAL 1.9dev, released 2011/01/18" LIBXML="2.7.8" USE_STATS"

>6) Your query should probably more look like (if you have a recent version of PostGIS):
> SELECT ST_MapAlgebraExpr(a.rast, 4, b.rast, 3, '([rast1] - [rast2]) / ([rast1] + [rast2])::float')
>FROM nclraster1 a, nclraster1 b;

Attempted and returns: "Neither raster provided has a NODATA value for the specified band indices.  NODATA value set to minimum possible for 8BUI"

ERROR:  syntax error at or near "["
LINE 1: SELECT (([$1] - [$2]) / ([$1] + [$2])::FLOAT)::double precis...
                 ^
QUERY:  SELECT (([$1] - [$2]) / ([$1] + [$2])::FLOAT)::double precision

Regards,
James
GIS Undergraduate
Reply | Threaded
Open this post in threaded view
|

Re: NDVI Calculation from two bands within one Raster

Pierre Racine-2
> > 5) What is the result of:SELECT PostGIS_Full_Version();
> "POSTGIS="2.0.0SVN" GEOS="3.3.2dev-CAPI-1.7.2" PROJ="Rel. 4.6.1, 21 August
> 2008" GDAL="GDAL 1.9dev, released 2011/01/18" LIBXML="2.7.8" USE_STATS"
>
> >6) Your query should probably more look like (if you have a recent
> >version
> of PostGIS):
> > SELECT ST_MapAlgebraExpr(a.rast, 4, b.rast, 3, '([rast1] - [rast2]) /
> >([rast1] + [rast2])::float') FROM nclraster1 a, nclraster1 b;
>
> Attempted and returns: "Neither raster provided has a NODATA value for the
> specified band indices.  NODATA value set to minimum possible for 8BUI"
>
> ERROR:  syntax error at or near "["
> LINE 1: SELECT (([$1] - [$2]) / ([$1] + [$2])::FLOAT)::double precis...
>                  ^
> QUERY:  SELECT (([$1] - [$2]) / ([$1] + [$2])::FLOAT)::double precision

You don't have a recent version (we are at beta 4), so try this instead (just remove the braquets):

SELECT ST_MapAlgebraExpr(a.rast, 4, b.rast, 3, '(rast1 - rast2) /(rast1 + rast2)::float') FROM nclraster1 a, nclraster1 b;

Pierre
_______________________________________________
postgis-users mailing list
[hidden email]
http://postgis.refractions.net/mailman/listinfo/postgis-users
Reply | Threaded
Open this post in threaded view
|

Re: NDVI Calculation from two bands within one Raster

JamesH
Attempted
SELECT ST_MapAlgebraExpr(a.rast, 4, b.rast, 3, '(rast1 - rast2) /(rast1 + rast2)::float') FROM nclraster1 a, nclraster1 b;

Returned "Neither raster provided has a NODATA value for the specified band indices.  NODATA value set to minimum possible for 8BUI"

Updated NoDataValues for bands 3 and 4 to 255 (were previously blank) and when re-ran, returns:
"ERROR:  division by zero
CONTEXT:  SQL statement "SELECT (($1 - $2) /($1 + $2)::FLOAT)::double precision"

Regards,
James
GIS Undergraduate
Reply | Threaded
Open this post in threaded view
|

Re: NDVI Calculation from two bands within one Raster

Pierre Racine-2
> Attempted
> SELECT ST_MapAlgebraExpr(a.rast, 4, b.rast, 3, '(rast1 - rast2) /(rast1 +
> rast2)::float') FROM nclraster1 a, nclraster1 b;
>
> Returned "Neither raster provided has a NODATA value for the specified band
> indices.  NODATA value set to minimum possible for 8BUI"
>
> Updated NoDataValues for bands 3 and 4 to 255 (were previously blank) and
> when re-ran, returns:
> "ERROR:  division by zero
> CONTEXT:  SQL statement "SELECT (($1 - $2) /($1 + $2)::FLOAT)::double
> precision"

Then you must have some 0s in band 3? What is the result of

SELECT ST_ValueCount(rast, 3, 0.0)
FROM nclraster1

What do you want/expect as result when band 3 = 0?
_______________________________________________
postgis-users mailing list
[hidden email]
http://postgis.refractions.net/mailman/listinfo/postgis-users
Reply | Threaded
Open this post in threaded view
|

Re: NDVI Calculation from two bands within one Raster

JamesH
Apologies, this was an grave error on my part as I didn't account for areas of the raster that have zero data values from where it has been georectified.

The result was that I have 16742 zero values in band 3 and 16627 in band 4 so I need to introduce something that excludes zero values from the NDVI calculation.

Is there a simple way of doing this as I seem to be over-complicating the calculation?

Regards,
James
GIS Undergraduate
Reply | Threaded
Open this post in threaded view
|

Re: NDVI Calculation from two bands within one Raster

Pierre Racine-2
> Apologies, this was an grave error on my part as I didn't account for areas of the
> raster that have zero data values from where it has been georectified.

They are zeros or nodata? If they are nodata setting them as nodata with ST_SetbandnodataValue() should be sufficient for ST_MapAlgebraExpr() to ignore them.

> The result was that I have 16742 zero values in band 3 and 16627 in band 4 so I
> need to introduce something that excludes zero values from the NDVI
> calculation.
>
> Is there a simple way of doing this as I seem to be over-complicating the
> calculation?

As you might have understood, the problem does not really comes from the fact that some values are equal to 0, but that the sum is equal to 0...

If they are zeros (not nodata), just make a CASE in the expression:

'WHEN (rast1 + rast2 = 0) THEN XXX ELSE (rast1 - rast2) /(rast1 + rast2)::float END'

Pierre

_______________________________________________
postgis-users mailing list
[hidden email]
http://postgis.refractions.net/mailman/listinfo/postgis-users
Reply | Threaded
Open this post in threaded view
|

Re: NDVI Calculation from two bands within one Raster

Bryce L Nordgren


'WHEN (rast1 + rast2 = 0) THEN XXX ELSE (rast1 - rast2) /(rast1 + rast2)::float END'

> Returned "Neither raster provided has a NODATA value for the specified band
> indices.  NODATA value set to minimum possible for 8BUI"

Ummm,

NDVI is calculated on reflectance values in these bands, not raw DN. Your inputs should already be floating point values, and they should never be less than 0.0 or greater than 1.0. The big concern here is that each band is likely to have a different relationship between DN and reflectance, making the above equation produce something...unexpected.




This electronic message contains information generated by the USDA solely for the intended recipients. Any unauthorized interception of this message or the use or disclosure of the information it contains may violate the law and subject the violator to civil or criminal penalties. If you believe you have received this message in error, please notify the sender and delete the email immediately.

_______________________________________________
postgis-users mailing list
[hidden email]
http://postgis.refractions.net/mailman/listinfo/postgis-users
Reply | Threaded
Open this post in threaded view
|

Re: NDVI Calculation from two bands within one Raster

JamesH
This post was updated on .
'NDVI is calculated on reflectance values in these bands, not raw DN. Your inputs should already be floating point values, and they should never be less than 0.0 or greater than 1.0.'

Ok, my raw data was 8BUI which I went back and changed to 64 bit float to try get past this but still, using the equation above is returning:

"NOTICE:  Neither raster provided has a NODATA value for the specified band indices.  NODATA value set to minimum possible for 64BF
ERROR:  division by zero
CONTEXT:  SQL statement "SELECT (($1 - $2) /($1 + $2))::double precision"

I went back and re-edited the raw raster image to remove the black areas from where it has been geo-rectified but still can't seem to get to grips with what is going wrong in the calculation.

This is the query as it stands:
SELECT
        ST_MapAlgebraExpr(a.rast, 4, b.rast, 3, '(rast1 - rast2) /(rast1 + rast2)::float')
CASE
        WHEN (rast1 + rast2 = 0) THEN XXX
        ELSE (rast1 - rast2) /(rast1 + rast2)::float
END
FROM
        ncl2img a, ncl2img b;

Returns Error: syntax error at or near "CASE" - unsure as to why?
note: ncl2img is the edited 64BF version of the previous nclraster1
GIS Undergraduate
Reply | Threaded
Open this post in threaded view
|

Re: NDVI Calculation from two bands within one Raster

Pierre Racine-2
You did not really answer my previous questions...

You can easily identify the faulty pixels by changing the expression:

SELECT ST_ValueCount(ST_MapAlgebraExpr(a.rast, 4, b.rast, 3, 'rast1 +
rast2'), 1, 0.0) FROM nclraster1 a, nclraster1 b;

Pierre

> -----Original Message-----
> From: [hidden email] [mailto:postgis-users-
> [hidden email]] On Behalf Of JamesH
> Sent: Tuesday, March 27, 2012 7:28 AM
> To: [hidden email]
> Subject: Re: [postgis-users] NDVI Calculation from two bands within one Raster
>
> NDVI is calculated on reflectance values in these bands, not raw DN. Your inputs
> should already be floating point values, and they should never be less than 0.0 or
> greater than 1.0.
>
> Ok, my raw data was 8BUI which I went back and changed to 64 bit float to try
> get past this but still, using the equation above is returning:
>
> "NOTICE:  Neither raster provided has a NODATA value for the specified band
> indices.  NODATA value set to minimum possible for 64BF
> ERROR:  division by zero
> CONTEXT:  SQL statement "SELECT (($1 - $2) /($1 + $2))::double precision"
>
> I went back and re-edited the raw raster image to remove the black areas from
> where it has been geo-rectified but still can't seem to get to grips with what is
> going wrong in the calculation.
>
> --
> View this message in context: http://postgis.17.n6.nabble.com/NDVI-
> Calculation-from-two-bands-within-one-Raster-tp4656995p4660390.html
> Sent from the PostGIS - User mailing list archive at Nabble.com.
> _______________________________________________
> postgis-users mailing list
> [hidden email]
> http://postgis.refractions.net/mailman/listinfo/postgis-users
_______________________________________________
postgis-users mailing list
[hidden email]
http://postgis.refractions.net/mailman/listinfo/postgis-users
Reply | Threaded
Open this post in threaded view
|

Re: NDVI Calculation from two bands within one Raster

JamesH
Apologies Pierre, the black areas around the original image (below) were the nodata values from where it had been oreintated.


SELECT ST_ValueCount(rast, 3, 0.0) FROM nclraster1
returned 16742 and 16627 in bands 3 and 4 respectively.
and
SELECT ST_ValueCount(ST_MapAlgebraExpr(a.rast, 4, b.rast, 3, 'rast1 + rast2'), 1, 0.0) FROM nclraster1 a, nclraster1 b;
returned 16223 pixels.

Both bands 3 and 4 now have nodata values of 255.

How would you remove the nodatavalues from the calculation?

Kind Regards,
James
GIS Undergraduate
Reply | Threaded
Open this post in threaded view
|

Re: NDVI Calculation from two bands within one Raster

Pierre Racine-2
> SELECT ST_ValueCount(ST_MapAlgebraExpr(a.rast, 4, b.rast, 3, 'rast1 + rast2'),
> 1, 0.0) FROM nclraster1 a, nclraster1 b; returned 16223 pixels.
>
> Both bands 3 and 4 now have nodata values of 255.
>
> How would you remove the nodatavalues from the calculation?

Nodata values, as far as they are properly set in the raster, are NOT part of the computation.

1) Make sure the nodata values are properly set with ST_BandNodataValue()

2) You have to get rid or decide what to do with all those pixels where the sum is equal to 0.

Pierre
_______________________________________________
postgis-users mailing list
[hidden email]
http://postgis.refractions.net/mailman/listinfo/postgis-users
Reply | Threaded
Open this post in threaded view
|

Re: NDVI Calculation from two bands within one Raster

JamesH
Have upgraded to Beta4,
new error message returning for
'SELECT ST_MapAlgebraExpr(a.rast, 4, b.rast, 3, '(rast1 - rast2) /(rast1 + rast2)::float') FROM nclraster1 a, nclraster1 b;'

As  column "rast1" does not exist.
With the previous PostGIS raster version, this has not returned before.
Any advice?

Regards,
James
GIS Undergraduate
Reply | Threaded
Open this post in threaded view
|

Re: NDVI Calculation from two bands within one Raster

JamesH
When altering the query to be:
SELECT
        ST_MapAlgebraExpr(a.rast, 4, b.rast, 3, '(a.rast - b.rast) /(a.rast + b.rast)::float')
FROM
        nclraster1 a, nclraster1 b;

instead of rast1 - rast2 etc.

Returns an error of
ERROR:  missing FROM-clause entry for table "a"
LINE 1: SELECT ((a.rast - b.rast) /(a.rast + b.rast)::float)::double...

Can anyone explain this?

Kind Regards,
James
GIS Undergraduate
Reply | Threaded
Open this post in threaded view
|

Re: NDVI Calculation from two bands within one Raster

Pierre Racine-2
> SELECT
> ST_MapAlgebraExpr(a.rast, 4, b.rast, 3, '(a.rast - b.rast) /(a.rast +
> b.rast)::float')
> FROM
> nclraster1 a, nclraster1 b;
>
> instead of rast1 - rast2 etc.
>
> Returns an error of
> ERROR:  missing FROM-clause entry for table "a"
> LINE 1: SELECT ((a.rast - b.rast) /(a.rast + b.rast)::float)::double...
>
> Can anyone explain this?

In the expression, you don't refer to the pixels of a.rast with "a.rast" but with "rast1". so your query should be:

SELECT ST_MapAlgebraExpr(a.rast, 4, b.rast, 3, '(rast1 - rast2) /(rast1 + rast2)::float') FROM nclraster1 a, nclraster1 b;

Pierre
_______________________________________________
postgis-users mailing list
[hidden email]
http://postgis.refractions.net/mailman/listinfo/postgis-users
Reply | Threaded
Open this post in threaded view
|

Re: NDVI Calculation from two bands within one Raster

JamesH
Sorry to keep reviving this thread, espcially with a basic question.

I am using:
SELECT ST_MapAlgebraExpr(a.rast, 4, b.rast, 3, '([rast1] - [rast2]) /([rast1] + [rast2])::float')
WHEN (rast1 + rast2 = 0) THEN = 999 ELSE ([rast1] - [rast2]) /([rast1] + [rast2])::float END'
FROM ndvi a, ndvi b;

Where the WHEN statement is being used to prevent the error division by zero.
However this returns: syntax error at or near "WHEN"

Can anyone advise on where I'm structuring this query wrongly.

Regards,
James
GIS Undergraduate
Reply | Threaded
Open this post in threaded view
|

Re: NDVI Calculation from two bands within one Raster

Pierre Racine-2
Your CASE syntax is wrong:

http://www.postgresql.org/docs/9.1/interactive/functions-conditional.html#FUNCTIONS-CASE

Should look like:

ST_MapAlgebraExpr(a.rast, 4, b.rast, 3, 'CASE WHEN ([rast1] + [rast2]) = 0 THEN 999 ELSE ([rast1] - [rast2]) /([rast1] + [rast2])::float END')
FROM ndvi a, ndvi b;

Pierre

> -----Original Message-----
> From: [hidden email] [mailto:postgis-users-
> [hidden email]] On Behalf Of JamesH
> Sent: Wednesday, April 11, 2012 10:56 AM
> To: [hidden email]
> Subject: Re: [postgis-users] NDVI Calculation from two bands within one Raster
>
> Sorry to keep reviving this thread, espcially with a basic question.
>
> I am using:
> SELECT ST_MapAlgebraExpr(a.rast, 4, b.rast, 3, '([rast1] - [rast2]) /([rast1] +
> [rast2])::float') WHEN (rast1 + rast2 = 0) THEN = 999 ELSE ([rast1] - [rast2])
> /([rast1] + [rast2])::float END'
> FROM ndvi a, ndvi b;
>
> Where the WHEN statement is being used to prevent the error division by zero.
> However this returns: syntax error at or near "WHEN"
>
> Can anyone advise on where I'm structuring this query wrongly.
>
> Regards,
> James
>
> -----
> GIS Undergraduate
> --
> View this message in context: http://postgis.17.n6.nabble.com/NDVI-
> Calculation-from-two-bands-within-one-Raster-tp4656995p4832397.html
> Sent from the PostGIS - User mailing list archive at Nabble.com.
> _______________________________________________
> postgis-users mailing list
> [hidden email]
> http://postgis.refractions.net/mailman/listinfo/postgis-users
_______________________________________________
postgis-users mailing list
[hidden email]
http://postgis.refractions.net/mailman/listinfo/postgis-users
Reply | Threaded
Open this post in threaded view
|

Re: NDVI Calculation from two bands within one Raster

JamesH
Thanks, that has worked and can visualise through OpenJump.

The Raster it produces (below), has only values of zero or one in its
attributes.
I'm assuming in the calculation its not keeping the values as Floats but its
converting it to integer? Must be rounding them to either zero or one.

Is there a way to solve this i.e. keep the pixels as floating integers?

http://postgis.17.n6.nabble.com/file/n4832776/ndvi.jpg 

Regards,
James

-----
GIS Undergraduate
--
View this message in context: http://postgis.17.n6.nabble.com/NDVI-Calculation-from-two-bands-within-one-Raster-tp4656995p4832776.html
Sent from the PostGIS - User mailing list archive at Nabble.com.
_______________________________________________
postgis-users mailing list
[hidden email]
http://postgis.refractions.net/mailman/listinfo/postgis-users
GIS Undergraduate
Reply | Threaded
Open this post in threaded view
|

Re: NDVI Calculation from two bands within one Raster

JamesH
In reply to this post by Pierre Racine-2
Thanks, that has worked and can visualise through OpenJump.

The Raster it produces (below), has only values of zero or one in its attributes.
 I'm assuming in the calculation its not keeping the values as Floats but its converting it to integer? Must be rounding them to either zero or one.
 
Is there a way to solve this i.e. keep the pixels as floating integers?


Kind Regards,
James
GIS Undergraduate
12