issue_comments: 457584997

Here is an old chunk of code I wrote awhile back to do this. Please note three things. There is the metadata attached to the file (I think it was through "tags"), metadata attached to the metadata "meta" variable, and some metadata that is attached on a per-band basis. It can be problematic when you assume that the info is global to the image and is embedded somehow (it took me weeks to figure some of this out).
Also note that I do per-band and image statistics... Also, I did not keep good enough notes and cannot remember where I got some of the hints and are just as likely to come from published examples that have been hacked to marginally work. Also, the .xml weirdness has to do in part with historic artifacts of our particular dataset that is over 3.5 petabytes, and cannot easily be updated, and is easier to hack in the code.

Hope this helps:

=====================

def to_tiff(data, fname, template=None, **kwargs): import numpy as np

 # check and promote the number of dimentio(1)ns for consistency
 nbands = data.ndim
 if 2 == nbands:
     # expand the array so that it is least 3D (ie stacks of

surfaces) import numpy as np data = np.expand_dims(data,axis=0) elif 3 != nbands: # nothing to do if it is already 3D print("Error: to_tiff can only currently deal with 2D and 3D data") return

 profile = {}
 tags = {}
 tmpl = None
 if template:
     tmpl = rasterio.open(template,'r')
     profile = tmpl.profile.copy()
     tags = tmpl.tags()

 # the metadata should be appended.  Cache here to
 # simplify variable replacement below.
 meta = {}
 if 'meta' in profile:
     meta.update(profile['meta'])
 if 'meta' in kwargs:
     meta.update(kwargs['meta'])

 # overwrite anything inheritied from the template with
 # user supplied args
 profile.update(kwargs)

 # overwrite bits that write the array as geotiff and
 # save the cached metadata
 profile['driver'] = 'GTiff'
 profile['count'] = data.shape[0]
 profile['width'] = data.shape[2]
 profile['height'] = data.shape[1]
 profile['meta'] = meta

 if 'dtype' not in profile:
     profile['dtype'] = type(data[0,0,0])

 # if you do not remove the previously associated .xml file,
 # then the tags and metadata can get corrupted.
 try:
     os.remove(fname)
     os.remove(fname+".xml")
 except:
     pass

 # now create and save the array to a file
 with rasterio.open(fname,'w',**profile) as out:
     for b in range(data.shape[0]):
         #print("\nprocessing band %d"%(b+1))
         out.write(data[b].astype(profile['dtype']), b+1)

         # caluclate the stats for each band
         # not sure what the proper name for per band stats is in

QGIS stats = { 'STATISTICS_MINIMUM': np.nanmin(data[b]), 'STATISTICS_MAXIMUM': np.nanmax(data[b]), 'STATISTICS_MEAN': np.nanmean(data[b]), 'STATISTICS_STDDEV': np.nanstd(data[b])} out.update_tags(b+1,**stats) #print(" stats= %s"%str(stats))

     # now calculate the stats across all the bands
     stats = {
             'STATISTICS_MINIMUM': np.nanmin(data),
             'STATISTICS_MAXIMUM': np.nanmax(data),
             'STATISTICS_MEAN': np.nanmean(data),
             'STATISTICS_STDDEV': np.nanstd(data)}

     out.update_tags(**tags)
     if 'tags' in kwargs:
         out.update_tags(**kwargs['tags'])

     out.update_tags(**stats)
     #print("\n  overall stats= %s\n"%str(stats))

 del tmpl

On Jan 23 2019 1:29 PM, Guillaume Eynard-Bontemps wrote:

Thanks @djhoese @ebo.

@ebo if you have some examples, that would be really cool!