2

I produce NDVI data from Landsat8 (UTM). End of the calculation, I want to export result EPSG:4326 so I used code that is below. I realized my new coordinates and qgis gdal warp result are different. lrx and lry are coordinates calculated wrong. Any idea?

enter image description here

g = gdal.Open(input_path)        
red = g.ReadAsArray()
s_srs = g.GetProjectionRef()        
osng = osr.SpatialReference ()
osng.SetFromUserInput ( s_srs )
wgs84 = osr.SpatialReference ()
wgs84.ImportFromEPSG ( 4326 )
tx = osr.CoordinateTransformation(osng,wgs84)
geo_t = g.GetGeoTransform ()        
x_size = g.RasterXSize # Raster xsize
y_size = g.RasterYSize
(ulx, uly, ulz ) = tx.TransformPoint( geo_t[0], geo_t[3])
x=geo_t[0] + geo_t[1]*x_size
y=geo_t[3] + geo_t[5]*y_size
(lrx, lry, lrz ) = tx.TransformPoint(x,y)

#gdal TransformPoint result
38.839209332307 38.528922768342554 0.0
41.406489023769865 36.39800679336268 0.0
#qgis gdalwarp result
#Extent 38.8392093323070071,36.3980038478590018 : 
41.4749801893845103,38.5290331355726110

2 Answers 2

1

I am not sure if you considered that if your warp operation is rotating the image then the original top-left pixel is not at the top-left corner of the envelope of the warped image.

enter image description here

A test with gdalwarp and image LC08_L1TP_033033_20180706_20180717_01_T1_B5.tif

Command

gdalwarp -t_srs epsg:4326 LC08_L1TP_033033_20180706_20180717_01_T1_B5.tif warped.tif -co tiled=yes -co compress=deflate

Corner coordinates of the original and warped

Corner Coordinates (original):
Upper Left  (  413385.000, 4423215.000) (106d 0'50.49"W, 39d57'16.71"N)
Lower Left  (  413385.000, 4187985.000) (105d59' 3.56"W, 37d50' 6.66"N)
Upper Right (  645015.000, 4423215.000) (103d18' 8.98"W, 39d56'47.93"N)
Lower Right (  645015.000, 4187985.000) (103d21' 7.93"W, 37d49'39.97"N)

Corner Coordinates (warped):
Upper Left  (-106.0140260,  39.9590749) (106d 0'50.49"W, 39d57'32.67"N)
Lower Left  (-106.0140260,  37.8276359) (106d 0'50.49"W, 37d49'39.49"N)
Upper Right (-103.3023551,  39.9590749) (103d18' 8.48"W, 39d57'32.67"N)
Lower Right (-103.3023551,  37.8276359) (103d18' 8.48"W, 37d49'39.49"N)

The images fit perfectly in QGIS if transformed on-the-fly into same CRS. Still the corner coordinates in lon-lat degrees differ. The reason is that the EPSG:32613 corner coordinates, when transformed into EPSG:4326, do not make a south-north oriented rectangle. But when the corner coordinates are computed for image that is warped to EPSG:4326 the lon-lat coordinates do make a south-north oriented rectangle.

Here is a screen capture from the lower right corner of the area covered by the Landsat image. The purple rectangle is native EPSG:32613 bounding box re-projected into EPSG:4326. The orange rectangle is the EPSG:4326 bounding from the image that is re-projected into EPSG:4326.

If you open the both images, native and warped into QGIS and toggle the visibility you will see that the net image area fits perfectly but the black nodata area has a slightly different shape.

enter image description here

8
  • if you match the coordinates, the problem is calculation of lrx and lry
    – Iroh
    Commented Feb 13, 2020 at 15:33
  • I do not understand what you mean. What I mean is that Extent 38.8392093323070071,36.3980038478590018 : 41.4749801893845103,38.5290331355726110 is the extent of the red rectangle while your calculation gives the location of the yellow vertex at the corner of the rotated green rectangle. It is correct result that they do not match.
    – user30184
    Commented Feb 13, 2020 at 15:44
  • I said, gdal.warp and my code give the min and max lat-lon value. Also if you check the code I used x_size and y_size to calculate the left lat -lon(lrx,lry) so I think theoretically both results should be matched.
    – Iroh
    Commented Feb 13, 2020 at 16:07
  • Give the corner coordinates and native CRS of your input image.
    – user30184
    Commented Feb 13, 2020 at 17:44
  • My input image is Landsat8. With link, you can download sample data quickly
    – Iroh
    Commented Feb 13, 2020 at 21:32
0

I solved the problem with memory driver and gdal Warp. I share my code and my resources in below.

def calculate_ndvi(self,red_path,nir_path,output_path):
        g = gdal.Open(red_path)        
        red = g.ReadAsArray().astype('float')
        g = gdal.Open(nir_path)        
        nir = g.ReadAsArray().astype('float')

        numpy.seterr(divide='ignore', invalid='ignore')
        ndvi = numpy.where(((nir+red)==0.),-32768,(nir-red)/(nir+red))
        ndvi_scale=ndvi*10000
        print(output_path)
        #raw image ul coordinate  e.g (348285.0, 30.0, 0.0, 4265115.0, 0.0, -30.0)  
        s_srs = g.GetProjectionRef()
        osng = osr.SpatialReference ()
        osng.SetFromUserInput ( s_srs )
        geo_t = g.GetGeoTransform ()
        x_size = g.RasterXSize # Raster xsize
        y_size = g.RasterYSize # Raster ysize
        mem_drv = gdal.GetDriverByName( 'MEM' )
        dest = mem_drv.Create('', x_size,y_size, 1, gdal.GDT_Int16)
        dest.SetGeoTransform( geo_t )
        dest.SetProjection ( osng.ExportToWkt() )
        dest.GetRasterBand(1).SetNoDataValue(-32768)
        dest.GetRasterBand(1).WriteArray(ndvi_scale)
        gdal.Warp(output_path, dest, format = 'GTiff', dstSRS = 'EPSG:4326 ')   
        dst_ds=None

source1

source2

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.