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