I have some test code that is using rsgislib to segment an image. From rsgislib's documentation I run:
from rsgislib.segmentation import segutils
inputImg = 'clipped_geotiff_image.tif'
outputClumps = 'clipped_jers1palsar_stack_clumps_elim_final.tif'
outputMeanImg = 'clipped_jers1palsar_stack_clumps_elim_final_mean.tif'
segutils.runShepherdSegmentation(inputImg, outputClumps, outputMeanImg, minPxls=100)
The outputMeanImg does not contain the coordinate system and other spatial information after trying to explicitly add it by:
rsgislib.imageutils.copySpatialAndProjFromImage(outputMeanImg, inputImg)
or
rsgislib.imageutils.assignProj(outputMeanImg, rsgislib.imageutils.getWKTProjFromImage(outputMeanImg), None)
For details, I have provided the output of gdalinfo below. Can anyone suggest how to set the corrdinate system and other projection information or to convert the subdataset BAND1/DATA to a conventional raster file?
=====
% gdalinfo clip_jers1palsar_stack_clumps_elim_final_mean.kea
Driver: HDF5/Hierarchical Data Format Release 5 Files: clip_jers1palsar_stack_clumps_elim_final_mean.kea Size is 512, 512 Coordinate System is `' Metadata: BAND1_DATA_BLOCK_SIZE=256d
BAND1_DATA_CLASS=IMAGE BAND1_DATA_IMAGE_VERSION=1.2
BAND1_NO_DATA_VAL_NO_DATA_DEFINED=
BAND1_OVERVIEWS_OVERVIEW1_BLOCK_SIZE=256d
BAND1_OVERVIEWS_OVERVIEW1_CLASS=IMAGE
BAND1_OVERVIEWS_OVERVIEW1_IMAGE_VERSION=1.2
BAND1_OVERVIEWS_OVERVIEW2_BLOCK_SIZE=187d
BAND1_OVERVIEWS_OVERVIEW2_CLASS=IMAGE
BAND1_OVERVIEWS_OVERVIEW2_IMAGE_VERSION=1.2
BAND1_OVERVIEWS_OVERVIEW3_BLOCK_SIZE=93d
BAND1_OVERVIEWS_OVERVIEW3_CLASS=IMAGE
BAND1_OVERVIEWS_OVERVIEW3_IMAGE_VERSION=1.2
BAND1_OVERVIEWS_OVERVIEW4_BLOCK_SIZE=46d
BAND1_OVERVIEWS_OVERVIEW4_CLASS=IMAGE
BAND1_OVERVIEWS_OVERVIEW4_IMAGE_VERSION=1.2 Subdatasets:
SUBDATASET_1_NAME=HDF5:"clip_jers1palsar_stack_clumps_elim_final_mean.kea"://BAND1/ATT/DATA/FLOAT SUBDATASET_1_DESC=[256x1] //BAND1/ATT/DATA/FLOAT (64-bit floating-point)
SUBDATASET_2_NAME=HDF5:"clip_jers1palsar_stack_clumps_elim_final_mean.kea"://BAND1/DATA SUBDATASET_2_DESC=[2500x1500] //BAND1/DATA (16-bit unsigned integer)
SUBDATASET_3_NAME=HDF5:"clip_jers1palsar_stack_clumps_elim_final_mean.kea"://BAND1/OVERVIEWS/OVERVIEW1 SUBDATASET_3_DESC=[625x375] //BAND1/OVERVIEWS/OVERVIEW1 (16-bit unsigned integer)
SUBDATASET_4_NAME=HDF5:"clip_jers1palsar_stack_clumps_elim_final_mean.kea"://BAND1/OVERVIEWS/OVERVIEW2 SUBDATASET_4_DESC=[312x187] //BAND1/OVERVIEWS/OVERVIEW2 (16-bit unsigned integer)
SUBDATASET_5_NAME=HDF5:"clip_jers1palsar_stack_clumps_elim_final_mean.kea"://BAND1/OVERVIEWS/OVERVIEW3 SUBDATASET_5_DESC=[156x93] //BAND1/OVERVIEWS/OVERVIEW3 (16-bit unsigned integer)
SUBDATASET_6_NAME=HDF5:"clip_jers1palsar_stack_clumps_elim_final_mean.kea"://BAND1/OVERVIEWS/OVERVIEW4 SUBDATASET_6_DESC=[78x46] //BAND1/OVERVIEWS/OVERVIEW4 (16-bit unsigned integer) Corner Coordinates: Upper Left ( 0.0, 0.0) Lower Left ( 0.0, 512.0) Upper Right ( 512.0, 0.0) Lower Right ( 512.0, 512.0) Center ( 256.0, 256.0)
I have found a partial solution, but not in rsgislib.
Using GDAL, I can open a subdataset:
from osgeo import gdal
# get the projection and geotransform from the original input image
ids = gdal.Open(inputImg, gdal.GA_ReadOnly)
gt = ids.GetGeoTransform()
wkt = ids.GetProjectionRef()
print gt
del ids
# now set the projection and geotransform for the specific subdatsert
infile = 'HDF5:"'+outputMeanImg+'"://BAND1/DATA'
ods = gdal.Open(infile, gdal.GA_ReadOnly)
ods.SetProjection(wkt)
ods.SetGeoTransform(gt)
del ods
The above code sets the projection and geotransform for the subdataset I am interested in, but please note that I am opening the output as ReadOnly. If I set it to update, it fails, and it actually updates the selected dataset even though it is opened in read only mode. I consider that a bug, but it appears to work by creating a new .aux.xml file.
There may be an easer/better way to do this, but this works for now.