2

I am using gdal_rasterize to clip mask by shapefile. However, the parameter "-te" does not work.

I use subprocess to call GDAL util:

import subprocess

# define paths to raster and vector
inraster = r'I:\HY_project\HY1B-JPG2\2011\deepwater\test\H1B_OPER_OCT_L1A_20110120T_011319707_geo.img'
inshape = r'I:\HY_project\HY1B-JPG2\2011\deepwater\test\H1B_OPER_OCT_L1A_20110120T_011319707_geo.shp'

# the extent of shapefile
driver = ogr.GetDriverByName("ESRI Shapefile")
dataSource = driver.Open(shp, 0)
layer = dataSource.GetLayer()
feature = layer[0]
geom = feature.GetGeometryRef()
extent = geom.GetEnvelope()
xmin = extent[0]
xmax = extent[1]
ymin = extent[2]
ymax = extent[3]

cmd = 'gdal_rasterize -burn -999 -te '+str(xmin)[0:7]+' '+str(ymin)[0:7]+' '+str(xmax)[0:7]+' '+\
     str(ymax)[0:7]+' '+inshape + ' ' + inraster

cmd = 'gdal_rasterize -burn -999 -i -te 130.827 24.8857 140.357 26.8057 '+ inshape + ' ' + inraster

# run command
subprocess.call(cmd, shell=True)

If the parameter "-te" is not set then the following command line works:

cmd = 'gdal_rasterize -burn -999 -i '+ inshape + ' ' + inraster
  • what does gdalinfo say for the image? may be projection is different or lat/lon not lon/lat – Ian Turton Sep 17 at 15:03
  • 2
    Does image H1B_OPER_OCT_L1A_20110120T_011319707_geo.img already exist? What do you mean with "parameter -te does not work"? In which way it does not work? – user30184 Sep 17 at 15:30
0

I would create a temporary layer and use gdal.RasterizeLayer toburn a value to the raster. Create a masked array to set the nodata value.

    geom = feature.GetGeometryRef()
    name = feature.GetField("name")
    env = geom.GetEnvelope()
    min_x, min_y, max_x, max_y = env[0], env[2], env[1], env[3]

    minX = round(min_x / cellsize, 0) * cellsize
    maxX = round(max_x / cellsize, 0) * cellsize 
    minY = round(min_y / cellsize, 0) * cellsize
    maxY = round(max_y / cellsize) * cellsize 

    ncol = int((maxX - minX) / cellsize)
    nrow= int((maxY - minY) / cellsize)

    tmp_ds = gdal.GetDriverByName('MEM').Create('', ncol, nrow, 1, gdal.GDT_Byte)
    tmp_ds.SetGeoTransform((
        minX, cellsize, 0,
        maxY, 0, -cellsize,
    ))

    #  Create for target raster the same projection as for the value raster
    dsmem = driver.CreateDataSource("/vsimem/temporarely.shp")
    l_out_mem = dsmem.CreateLayer("memlayer", srs, ogr.wkbPolygon)
    # Add an ID field
    idField = ogr.FieldDefn("id", ogr.OFTInteger)
    l_out_mem.CreateField(idField)

    # Create the feature and set values
    featureDefn = l_out_mem.GetLayerDefn()
    out_feature = ogr.Feature(featureDefn)
    out_feature.SetField("id", 1)
    out_feature.SetGeometry(geom.Clone())
    l_out_mem.CreateFeature(out_feature)

    # burn value 1 to new raster    
    tmp_ds.SetProjection(srs.ExportToWkt())
    gdal.RasterizeLayer(tmp_ds, [1], l_out_mem, burn_values=[1])

    dsmem = None
    # clip same extent from srcImage
    cellX = int((minX - xOrigin) / pixelWidth)
    cellY = int((yOrigin - maxY) / -pixelHeight)

    nparr = rb.ReadAsArray(cellX, cellY, ncol, nrow)


    # create mask from rasterized layer
    bandmask = tmp_ds.GetRasterBand(1)
    datamask = bandmask.ReadAsArray(0, 0, ncol, nrow).astype(np.float)
    datamask[datamask == 0] = 'nan'

    # Mask zone of raster
    if (nparr is not None):
        zonedata = nparr * datamask
        zonedata[zonedata == src_ds_nd] = 'nan'
        n = np.count_nonzero(~np.isnan(zonedata))

    if (nparr is None or n < 1) :
        QMessageBox.information(None, "No data is found in feature: "+name,"")
        continue

    # set NoData as nan
    zonedata = zonedata.astype(np.float)
    zonedata[zonedata == src_ds_nd] = 'nan'

    # write the clipped raster to file   
    target_ds = gdal.GetDriverByName('GTiff').Create(os.path.join(OUT_FOLDER, name+'.tif'), nkol, nrij, 1, gdal_datatype)
    target_ds.SetGeoTransform((minX, cellsize, 0, maxY, 0, -cellsize, ))

    # Create for target raster the same projection as for the value raster
    target_ds.SetProjection(srs.ExportToWkt())
    outBand = target_ds.GetRasterBand(1)

    # write the data
    outBand.WriteArray(zonedata, 0, 0)
    outBand.FlushCache()
    outBand.SetNoDataValue(-9999)
    target_ds = None
    outBand = None

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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