title: Crop satellite image image based on a historical image with OpenCV in Python
body: <p>I have the following problem, I have a pair of two images one historical and one present-day satellite image and as the historical image covers a smaller area I want to crop the satellite images. Here one can see the code I wrote for this:</p>
  
  <pre><code>import numpy as np
  import cv2
  import os
  import imutils
  import math
  
  entries = os.listdir('../')
  
  refImage = 0
  histImages = []
  def loadImage(index):
      referenceImage = cv2.imread("../" + 'ref_' + str(index) + '.png')
  
      top = int(0.5 * referenceImage.shape[0])  # shape[0] = rows
      bottom = top
      left = int(0.5 * referenceImage.shape[1])  # shape[1] = cols
      right = left
  
      referenceImage = cv2.copyMakeBorder(referenceImage, top, bottom, left, right, cv2.BORDER_CONSTANT, None, (0,0,0))
      counter = 0
      for entry in entries:
          if entry.startswith("image_"+str(index)):
              refImage = referenceImage.copy()
              histImage = cv2.imread("../" + entry)
              #histImages.append(img)
              points = np.loadtxt("H2OPM/"+"CP_"+ entry[6:9] + ".txt", delimiter=",")
  
              vector_image1 = [points[0][0] - points[1][0], points[0][1] - points[1][1]] #hist
              vector_image2 = [points[0][2] - points[1][2], points[0][3] - points[1][3]] #ref
  
              angle = angle_between(vector_image1, vector_image2)
  
              hhist, whist, chist = histImage.shape
  
              rotatedImage = imutils.rotate(refImage, angle)
  
              x = int(points[0][2] - points[0][0])
              y = int(points[1][2] - points[1][0])
  
              crop_img = rotatedImage[x+left:x+left+hhist, y+top:y+top+whist]
  
              print("NewImageWidth:", (y+top+whist)-(y+top),(x+left+hhist)-(x+left))
  
              print(entry)
              print(x,y)
  
  
              counter += 1
              #histImage = cv2.line(histImage, (points[0][0], ), end_point, color, thickness) 
              cv2.imwrite("../matchedImages/"+'image_' + str(index) + "_" + str(counter) + '.png'  ,histImage)
              #rotatedImage = cv2.line(rotatedImage, (), (), (0, 255, 0), 9) 
              cv2.imwrite("../matchedImages/"+'ref_' + str(index) + "_" + str(counter) + '.png'  ,crop_img)
  </code></pre>
  
  <p>First, I load the original satellite image and pad it so I don't lose information due to the rotation, second, I load one of the matched historical images as well as the matched keypoints of the two images (i.e. a list of x_hist, y_hist, x_present_day, y_present_day). Third, I compute the rotation angle between the two images (which works) and fourth, I crop the image (and fifth, I save the images).</p>
  
  <p><strong>Problem</strong>: As stated the rotation works fine, but my program ends up cropping the wrong part of the image.</p>
  
  <p>I think that, due to the rotation, the boundaries (i.e. left, right, top, bottom) are no longer correct and I think this is where my problem lies, but I am not sure how to fix this problem.</p>
  
  <p><strong>Information that might help:</strong></p>
  
  <ul>
  <li>The images are both scaled the same way (so one pixel = approx. 1m)</li>
  <li>I have at least 6 keypoints for each image</li>
  </ul>