;;
;; geotiff file name, lon latを入力して対応する 画像座標(x, y)を返す
;;
function geo_transform, ifname, lon, lat, dlon=dlon, dlat=dlat, tile_size_x=tile_size_x, tile_size_y=tile_size_y
if n_elements(tile_size_x) eq 0 then begin
tile_size_x = 1
endif
if n_elements(tile_size_y) eq 0 then begin
tile_size_y = 1
endif
if n_elements(dlon) eq 0 then begin
dlon = 0d
endif
if n_elements(dlat) eq 0 then begin
dlat = 0d
endif
lon_arr = dindgen_x(tile_size_x,tile_size_y)*dlon + lon
lat_arr = dindgen_x(tile_size_x,tile_size_y,/y)*dlat + lat
;; sample ;;
input_image = query_tiff(ifname,im_info, geotiff=geotag)
;; 画像サイズ ;;
s = im_info.dimensions
;; UTM absolute coordinate ;;
xscale = geotag.ModelPixelScaleTag[0]
yscale = geotag.ModelPixelScaleTag[1]
tp = geotag.ModelTiePointTag[3:4]
;; tie point = (0, ny)
xOrigin = tp[0]
yOrigin = tp[1] ;- (yscale * s[1])
;; 必要なら 0.5pixずらす
xOrigin = xOrigin + (xscale * 0.5)
yOrigin = yOrigin + (yscale * 0.5)
mapCoord = cgGeoMap(ifname)
mapStruct = mapCoord -> GetMapStruct()
;; Lon, Lat => UTM x,y
utm_xy = MAP_PROJ_FORWARD(lon_arr, lat_arr, MAP_STRUCTURE=mapStruct)
;; UTM x,y => Image x, y (1 dimension) ;;
;; tie point = (0, ny)
tmp_im_x = (utm_xy[0,*] - xOrigin)/xscale
tmp_im_y = (yOrigin - utm_xy[1,*])/xscale
;; 1 dimension -> 2 dimension ;;
im_x = Reform(tmp_im_x[*], tile_size_x, tile_size_y)
im_y = Reform(tmp_im_y[*], tile_size_x, tile_size_y)
im_xy = dblarr(2,tile_size_x, tile_size_y)
im_xy[0,*,*] = im_x
im_xy[1,*,*] = im_y
return,im_xy
end