mask and extract cell values from a vrt file?

Question

I have raster data for built up areas around the globe with 40m resolution as vrt file, download data from here , and I am trying to crop the data by a mask and then extract color index value for each cell.

Note: another 2 files exist with the data: vrt.clr and vrt.ovr

Here is a sample of data: view of vrt data in arcmap.

My question: why I am getting empty cells values when I crop by mask ?

I have tried the following:

• extract by mask using arcmap toolbox

• using gdal in python 2.7

  import gdal ds = gdal.Open('input.vrt') ds = gdal.Translate('output.vrt', ds, projWin = [80.439,5.341,81.048,4.686]) ds = None

• I have also try to save the data as tif

Also, is there any way to read the color index value at given coordinates (x,y) after masking the data?

Answer 1

The data appears to be in the Pseudo Mercator projection (EPSG 3857). So therefore you should either specify the extent for projWin in that coordinate system, or add projWinSRS if you want to provide them in a different coordinate system.

Also, if you want gdal.Translate to output to a VRT file, you should add format='VRT. Because in your code snippet outputs to the default file format, which is GeoTIFF.

When i assume your coordinates are WGS84 (EPSG 4326), it defines a small region over the ocean south of Sri Lanka. That doesn't make much sense given the nature of the data.

If you want to read the array given by your coordinates you could use:

invrt = 'GHS_BUILT_LDSMT_GLOBE_R2015B_3857_38_v1_0.vrt'
outfile = '/vsimem/tmpfile'

ds = gdal.Translate(outfile, invrt, projWin=[80.439, 5.341, 81.048, 4.686], projWinSRS='EPSG:4326')
data = ds.ReadAsArray()
ds = None
gdal.Unlink(outfile)

The plotted array looks like: