677def wkblify_raster_header(options, ds, level, ulp, xsize = None, ysize = None):
678 """Writes WKT Raster header based on given GDAL into HEX-encoded WKB."""
679 assert ds is not None, "Error: Missing GDAL dataset"
680 assert level >= 1
681 assert len(ulp) == 2, "Error: invalid upper-left corner"
682
683 if xsize is None or ysize is None:
684 assert xsize is None and ysize is None
685 xsize = ds.RasterXSize
686 ysize = ds.RasterYSize
687
688
689 gt = get_gdal_geotransform(ds)
690 ul = calculate_geoxy(gt, (ulp[0], ulp[1]))
691 rt_ip = ( ul[0], ul[1] )
692 rt_skew = ( gt[2], gt[4] )
693 rt_scale = ( gt[1] * level, gt[5] * level )
694
695
696
697
698
699
700 hexwkb = ''
701
702 hexwkb += wkblify('B', options.endian)
703
704 hexwkb += wkblify('H', options.version)
705
706 if options.band is not None and options.band > 0:
707 hexwkb += wkblify('H', 1)
708 else:
709 hexwkb += wkblify('H', ds.RasterCount)
710 check_hex(hexwkb, 5)
711
712 hexwkb += wkblify('d', rt_scale[0])
713 hexwkb += wkblify('d', rt_scale[1])
714 hexwkb += wkblify('d', rt_ip[0])
715 hexwkb += wkblify('d', rt_ip[1])
716 hexwkb += wkblify('d', rt_skew[0])
717 hexwkb += wkblify('d', rt_skew[1])
718 hexwkb += wkblify('i', options.srid)
719 check_hex(hexwkb, 57)
720
721 hexwkb += wkblify('H', xsize)
722 hexwkb += wkblify('H', ysize)
723 check_hex(hexwkb, 61)
724
725 logit("MSG: Georeference: px = %s -> ul = %s \tresolution = %s \trotation = %s\n" \
727 return hexwkb
728