0

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:

enter image description here

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()
  • Are your images rotated correctly? Images are Row 0, Col 0 at top left and continue downward... YCellSize is usually negative... Tiff files support a rotation parameter which can be a bit of a pain - but I haven't seen that set in many years though. – Michael Stimson Mar 26 '15 at 21:20
  • Definitely something is wrong with the rotation - when I run this on images with 0 rotation it works as it should. I am taking the rotation angle given for creating the frame which is in degrees, converting it to radians and then applying some math to it to calculate the world file values you can see the math in the script – kflaw Mar 26 '15 at 21:41
1

I found that I was inputting the dimensions and rotation value incorrectly and therefore not grabbing the correct coordinates for the upper left corner of the image. Also, though I don't understand yet why, the sign on the rotation value needs to be reversed for the world file. Here is the updated 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 write world files for: ')
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()

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.