I have centroids for map images which I would like to georeference quickly and rougly via python script, using arcpy. I know the dimensions of the images, the dpi, scale, rotation and the coordinate system to put them in. My script builds frames or rectangle extents based off the centerpoint, and then I am trying to build the world file by inputting the resolution and factoring in the rotation angle - i used this post as a guide - http://www.omg.unb.ca/~jonnyb/processing/geotiff_tifw_format.html. I am using an attribute value in the points layer to build out the world file names to match the image names. The script runs without error, but when I load my frames in and my images with world files into a map, the images and frames do not match up. The top right corner of the image lands where the top left corner of the frame is. Here's an image of what it looks like:
I also found this post How to generate a world file with rotation but I couldn't find an answer from it. What am I doing wrong? Here is my code:
import arcpy, os
w = float(raw_input('provide approximate width of map frame in desired map units: '))/2
l = float(raw_input('provide approximate length of map frame in desired map units: '))/2
inputAngle = int(raw_input('provide angle for frame orientation: '))
pointLayer = r'C:\Temp\temp.gdb\QD_Georef_test_geocoded_CCS27z3e'
pointFields = ['POINT_X','POINT_Y']
polygonFCpath = r'C:\Temp\temp.gdb'
polygonFCname = 'test'
polygonFC = os.path.join(polygonFCpath,polygonFCname)
features = []
polyArray = {}
arcpy.env.overwriteOutput = True
arcpy.AddXY_management(pointLayer)
frameID = str(raw_input('Enter field name containing Frame ID to match to image names: '))
pointFields.insert(0,frameID)
polygonFields = [frameID,'UpperLeftX', 'UpperLeftY','SHAPE@']
srID = raw_input('Enter desired spatial reference for frames by code or name: ')
if srID[0].isdigit():
srInput = int(srID)
else:
srInput = srID
sr = arcpy.SpatialReference(srInput)
arcpy.CreateFeatureclass_management(polygonFCpath,polygonFCname,'POLYGON')
arcpy.DefineProjection_management(polygonFC,sr)
arcpy.AddField_management(polygonFC,frameID,'TEXT','','',50)
arcpy.AddField_management(polygonFC,'UpperLeftX','DOUBLE')
arcpy.AddField_management(polygonFC,'UpperLeftY','DOUBLE')
def RotateXY(x,y,xc=0,yc=0,angle=0,units='DEGREES'):
import math
x = float(x) - xc
y = float(y) - yc
if units == 'DEGREES':
angle = math.radians(angle)
xr = (x * math.cos(angle)) - (y * math.sin(angle)) + xc
yr = (x * math.sin(angle)) + (y * math.cos(angle)) + yc
return xr,yr
for row in arcpy.da.SearchCursor(pointLayer,pointFields):
featureList = []
upperLeftX = str(row[1] - w)
upperLeftY = str(row[2] + l)
upperRightX = str(row[1] + w)
upperRightY = str(row[2] + l)
lowerLeftX = str(row[1] - w)
lowerLeftY = str(row[2] - l)
lowerRightX = str(row[1] + w)
lowerRightY = str(row[2] - l)
upperLeftVal = RotateXY(upperLeftX,upperLeftY,row[1],row[2],inputAngle)
upperRightVal = RotateXY(upperRightX,upperRightY,row[1],row[2],inputAngle)
lowerLeftVal = RotateXY(lowerLeftX,lowerLeftY,row[1],row[2],inputAngle)
lowerRightVal = RotateXY(lowerRightX, lowerRightY,row[1],row[2],inputAngle)
shtID = row[0]
coordsList = []
coordsList.append(upperLeftVal)
coordsList.append(upperRightVal)
coordsList.append(lowerRightVal)
coordsList.append(lowerLeftVal)
polyArray[shtID, upperLeftVal[0],upperLeftVal[1]] = coordsList
for key, values in polyArray.iteritems():
features.append((key[0],key[1],key[2],arcpy.Polygon(arcpy.Array([arcpy.Point(*value) for value in values]),
arcpy.SpatialReference(srInput))))
with arcpy.da.InsertCursor(polygonFC,polygonFields) as cur:
for feature in features:
cur.insertRow(feature)
cellSize = float(raw_input('Provide cell size for georeferenced image (scale/dpi): '))
angleRadians = math.radians(inputAngle)
worldFileLine1 = cellSize * math.cos(angleRadians)
worldFileLine2 = -cellSize * math.sin(angleRadians)
worldFileLine3 = -cellSize * math.sin(angleRadians)
worldFileLine4 = -cellSize * math.cos(angleRadians)
imageWkspc = raw_input('Provide path location for images to be georeferenced: ')
arcpy.env.workspace = imageWkspc
imgList = arcpy.ListRasters()
for i in imgList:
for row in arcpy.da.SearchCursor(polygonFC,polygonFields):
mapSet = row[0].split('-')[0].lower() +'_'+row[0].split('-')[1].split('.')[0].zfill(5)+'_'+ row[0].split('.')[1].zfill(3)
if mapSet in i:
iName = i.split('.')[0]
worldFile = open(os.path.join(imageWkspc,iName+'.tfw'),'w')
worldFileLine5 = row[1]
worldFileLine6 = row[2]
print 'writing out world file for',i,'...'
worldFile.write(str(worldFileLine1)+'\n'+str(worldFileLine2)+'\n'+str(worldFileLine3)+'\n'+str(worldFileLine4)+'\n'+str(worldFileLine5)+'\n'+str(worldFileLine6))
worldFile.close()
