Community
Select to view content in your preferred language

Rotated raster cell size issue, ArcGIS does not read the metadata? (ASTER geotiff)

311
0
03-05-2026 07:16 PM
Labels (1)
Labels
SOP147
by
New Contributor
‎03-05-2026 07:16 PM

Hello. Thank you in advance to anyone helping with this matter.

 

Basically I was looking at some VNIR ASTER data (ASTER L2 reflectance so AST_07). I often use ENVI to take a first look at the data. On ENVI it looks good; it says projection: UTM 33N, 15m pixel size, rotation 10°. All pixels are aligned N-S. This is perfect
But when opening the data on ArcGIS Pro, the software does recognize the right geographic and projected system UTM 33N, but it doesn't handle the rotation in the metadata at all. As a result, I get a pixel size of 16.97 x 17.87 m, and rotation was not applied so that pixels are displayed slanted.
My issue is that I want to avoid resampling, because the value of each pixel is important and any resampling would change that.

Here is the detail of what is going on:
I used python to get me the file metadata without any influence from either ENVI nor ArcGIS
Most relevant info are:
• Top-left X: 464956.94075709174
• Pixel width (X): 14.769436715652846
• Rotation X: -2.6198738714536747
• Top-left Y: 4299748.463908005
• Rotation Y: -2.6198738714536747
• Pixel height (Y): -14.769436715652846
• Projection: PROJCS["WGS 84 / UTM zone 33N",GEOGCS["WGS 84" ............
• Columns: 4980
• Rows: 4200

ENVI knows to read that properly so that pixel size is calculated using the "pixel size along the pixel axes", and so true size = √(14.77² + 2.62²) = 15m ✅ (14.7 m to the East, 2.62 m to the South)
Similarly, rotation angle is θ = arctan(14.769/2.6199) = 10.06°, this is exactly what ENVI is giving me,

Now, in ArcGIS, the "metadata" gives:
• Cell Size X: 16.979m
• Cell Size Y: 17.876m
Well, coincidentally...
• X_extent = (Columns × Pixel width) + (Rows × Rotation) = (4980×14.769)+(4200×2.620) = 84,553.62 m
• Y_extent = (Rows × Pixel width) + (Columns × Rotation) = (4200×14.769) + (4980×2.620) = 75,077.4 m
Which resolves as:
• Cell Size X = X_extent / Columns = 16.978
• Cell Size Y = Y_extent / Rows = 17.876
ArcGIS is using the whole bounding box including the NaN edges created by the rotation to compute the pixel size. It is just a fundamentally different way to read the metadata.

I have no idea how to handle that. Is there a way in ArcGIS to make it read the metadata correctly with the rotation coefficients?

Here is a code to help with visualization of the rotation issue

import matplotlib.pyplot as plt
import matplotlib.patches as patches
import numpy as np

fig, ax = plt.subplots(figsize=(8, 7))
ax.set_xlim(0, 560)
ax.set_ylim(0, 480)
ax.set_aspect('equal')
ax.axis('off')
ax.invert_yaxis()  # SVG-style: y increases downward

# --- Rotated image corners ---
# angle = 10.06°, cos=0.9846, sin=0.1746
# Origin top-left: ox=100, oy=80, imgW=230, imgH=190
ox, oy = 100, 80
imgW, imgH = 230, 190
angle = np.radians(10.06)
c, s = np.cos(angle), np.sin(angle)

TL = np.array([ox, oy])
TR = np.array([ox + imgW*c, oy + imgW*s])
BR = np.array([ox + imgW*c - imgH*s, oy + imgW*s + imgH*c])
BL = np.array([ox - imgH*s, oy + imgH*c])

corners = np.array([TL, TR, BR, BL, TL])  # close polygon

# Bounding box
minX = min(TL[0], TR[0], BR[0], BL[0])
maxX = max(TL[0], TR[0], BR[0], BL[0])
minY = min(TL[1], TR[1], BR[1], BL[1])
maxY = max(TL[1], TR[1], BR[1], BL[1])

# --- Draw bounding box ---
bbox = patches.Rectangle((minX, minY), maxX-minX, maxY-minY,
                          linewidth=2, edgecolor='#ef4444',
                          facecolor='#fee2e2', linestyle='--', zorder=1)
ax.add_patch(bbox)

# --- Draw rotated image ---
img_poly = patches.Polygon(corners[:-1], closed=True,
                            facecolor='#bfdbfe', edgecolor='#1d4ed8',
                            linewidth=2, alpha=0.9, zorder=2)
ax.add_patch(img_poly)

# --- Grid lines on rotated image ---
for i in [1, 2, 3, 4]:
    t = i / 5
    p1 = TL + t * (TR - TL)
    p2 = BL + t * (BR - BL)
    ax.plot([p1[0], p2[0]], [p1[1], p2[1]], color='#1d4ed8', lw=0.6, alpha=0.4, zorder=3)
for i in [1, 2, 3]:
    t = i / 4
    p1 = TL + t * (BL - TL)
    p2 = TR + t * (BR - TR)
    ax.plot([p1[0], p2[0]], [p1[1], p2[1]], color='#1d4ed8', lw=0.6, alpha=0.4, zorder=3)

# --- North arrow ---
ax.annotate('', xy=(33, 30), xytext=(33, 55),
            arrowprops=dict(arrowstyle='->', color='#333', lw=1.8))
ax.text(29, 70, 'N', fontsize=13, color='#333')

# --- Rotation angle label ---
ax.text(133, 76, 'θ=10°', fontsize=9, color='#ea580c')

# --- One pixel decomposition at TL ---
# True 15m pixel vector (scaled: 1m ≈ 3.07px, so 15m ≈ 46px)
scale = 46  # visual pixel length
px = TL[0] + scale * c
py = TL[1] + scale * s

# True pixel (blue arrow)
ax.annotate('', xy=(px, py), xytext=(TL[0], TL[1]),
            arrowprops=dict(arrowstyle='->', color='#1d4ed8', lw=2.5), zorder=5)
ax.text(TL[0]+3, TL[1]-9, '15m (true pixel)', fontsize=9,
        color='#1d4ed8', fontweight='bold')

# East component gt1 (green dashed)
ax.annotate('', xy=(px, TL[1]), xytext=(TL[0], TL[1]),
            arrowprops=dict(arrowstyle='->', color='#16a34a', lw=1.8,
                           linestyle='dashed'), zorder=5)
ax.text((TL[0]+px)/2, TL[1]-6, 'gt₁=14.77m',
        fontsize=8, color='#16a34a', ha='center')

# South component gt2 (orange dashed)
ax.annotate('', xy=(px, py), xytext=(px, TL[1]),
            arrowprops=dict(arrowstyle='->', color='#ea580c', lw=1.8,
                           linestyle='dashed'), zorder=5)
ax.text(px+4, (TL[1]+py)/2, '|gt₂|=2.62m',
        fontsize=8, color='#ea580c')

# Right angle marker
sq = 6
ax.add_patch(patches.Rectangle((px-sq, TL[1]), sq, sq,
             fill=False, edgecolor='#888', lw=1, zorder=5))

# --- Bounding box width annotation ---
y_ann = maxY + 22
ax.annotate('', xy=(maxX, y_ann), xytext=(minX, y_ann),
            arrowprops=dict(arrowstyle='<->', color='#ef4444', lw=1.5))
ax.text((minX+maxX)/2, y_ann+14, 'X_extent = cols×gt₁ + rows×|gt₂|',
        fontsize=9, color='#ef4444', ha='center')

# --- Bounding box height annotation ---
x_ann = maxX + 20
ax.annotate('', xy=(x_ann, maxY), xytext=(x_ann, minY),
            arrowprops=dict(arrowstyle='<->', color='#ef4444', lw=1.5))
ax.text(x_ann+4, (minY+maxY)/2 - 8,  'Y_extent =', fontsize=9, color='#ef4444')
ax.text(x_ann+4, (minY+maxY)/2 + 5,  'rows×gt₁', fontsize=9, color='#ef4444')
ax.text(x_ann+4, (minY+maxY)/2 + 18, '+ cols×|gt₂|', fontsize=9, color='#ef4444')

# --- Labels ---
ax.text(minX+2, minY+14, 'Bounding box (ArcGIS reads this)',
        fontsize=9, color='#ef4444', fontweight='bold')
ax.text(ox+50, oy+110, 'Rotated image\n(true data)',
        fontsize=9, color='#1d4ed8', fontweight='bold', ha='center')

# Empty corners
ax.text(minX+2, maxY-10, 'empty corner', fontsize=8, color='#ef4444', alpha=0.7)
ax.text(TR[0]-10, minY+12, 'empty corner', fontsize=8, color='#ef4444', alpha=0.7)

# --- Formulas at bottom ---
ax.text(80, 410, '√(14.77² + 2.62²) = 15.0m  ← true pixel size (ENVI)',
        fontsize=10, color='#1d4ed8')
ax.text(80, 430, 'X_extent ÷ cols = 16.98m      ← what ArcGIS reports',
        fontsize=10, color='#ef4444')

plt.tight_layout()
plt.savefig('bounding_box.png', dpi=150, bbox_inches='tight')
plt.show()
0 Kudos
0 Replies