2

I'm trying to mosaic several rasters with gdal.BuildVRT in Python 2.7 (Windows system), but code runs with logical problems.

When I try to merge 3 Landsat images datasets, resulting image is a tile (Float32) that include 'nan' values in raster areas of two datasets (see next image).

Landsat Tile Float 32

When I do a mosaic with GdalBuildVRT in Python, it shows me the following problem:

Fail mosaic

The code used is the following:

from osgeo import gdal
import osr, numpy as np

def save_raster ( output_name, dataset, driver ,NaN_Value):
    """
    A function to save a 1-band raster using GDAL to the file indicated
    by ``output_name``. It requires a GDAL-accesible dataset to collect 
    the projection and geotransform.
    """
    # Open the reference dataset
    g = ( dataset )
    # Get the Geotransform vector
    geo_transform = g.GetGeoTransform ()
    x_size = g.RasterXSize # Raster xsize
    y_size = g.RasterYSize # Raster ysize
    srs = g.GetProjectionRef () # Projection
    raster_data = g.ReadAsArray()
    NaN_rast = NaN_Value
    # raster_data[raster_data == NaN_rast] = 'NaN'
    raster_data[raster_data == NaN_rast] = np.NaN
    # Need a driver object. By default, we use GeoTIFF
    driver = gdal.GetDriverByName ( driver )
    dataset_out = driver.Create ( output_name, x_size, y_size, 1, \
            gdal.GDT_Float32 )
    dataset_out.SetGeoTransform ( geo_transform )
    dataset_out.SetProjection ( srs )
    dataset_out.GetRasterBand ( 1 ).WriteArray ( \
            raster_data.astype(np.float32) )

TILE1 = "C:\Users\italo\Documents\GeoPRADS\LANDSAT\NUEVAS IMGS LANDSAT\Input ETa\TILE1.TIF" 
TILE2 = "C:\Users\italo\Documents\GeoPRADS\LANDSAT\NUEVAS IMGS LANDSAT\Input ETa\TILE2.TIF"
TILE3 = "C:\Users\italo\Documents\GeoPRADS\LANDSAT\NUEVAS IMGS LANDSAT\Input ETa\TILE3.TIF"    

in_subruta = TILE1+'-'+TILE2+'-'+TILE3 #string 
subruta = in_subruta.split('-') #to list
gdal.BuildVRTOptions(VRTNodata= 0)
b1 = gdal.BuildVRT("dummy", subruta,VRTNodata= "NaN" )        
save_raster(pathout+'/mosaic.tif',b1,"GTiff",0)

VRTNodata take the 'nan' data as a value.

So, where is the problem in the code?

3
  • sorry, the object pathout is a string value pathout = "C:\Users\italo\Documents\GeoPRADS\LANDSAT\NUEVAS IMGS LANDSAT\Input ETa" Commented Mar 17, 2017 at 19:28
  • Have you tried either b1 = gdal.BuildVRT("dummy", subruta, VRTNodata=np.NaN) or b1 = gdal.BuildVRT("dummy", subruta, VRTNodata=np.NaN, srcNodata=np.NaN)? Commented Mar 17, 2017 at 19:55
  • yeah, i tried but doesn't work Commented Mar 17, 2017 at 20:22

1 Answer 1

3

Ok, i solved the error

the value in srcNodata must be a string "nan"

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.