1

I am trying to make a 2D planer map of Mars, I have multiple non_georef images along with their geocentric lat and long. I am using gdal.GCP to georef these images; its working but something seems off as the mountains of Mars are not aligned in a way as they should be and all the ice is on the south pole and none on the north, is there something I am doing wrong?

Here is the code for reference

def get_georef(lat_path, long_path, image, image_out_geo):
    lat_mat = sio.loadmat(lat_path)
    long_mat = sio.loadmat(long_path)
lat = np.array(lat_mat['latitude'])
long = np.array(long_mat['longitude'])

dataset = gdal.Open(image)


shutil.copy(image, image_out_geo)

ds = gdal.Open(image_out_geo, gdal.GA_Update)

x1 = long[0][0]
y1 = lat[767][1023]

x2 = long[767][1023]
y2 = lat[0][0]

xm1 = long[0][1023]
ym1 = lat[767][0]

xm2 = long[767][0]
ym2 = lat[0][1023]

xc = long[383][511]
yc = lat[383][511]

xc1 = long[0][511]
yc1 = lat[0][511]

xc2 = long[383][0]
yc2 = lat[383][0]

xcm1 = long[383][1023]
ycm1 = lat[383][1023]

xcm2 = long[767][511]
ycm2 = lat[767][511]

gcps = [gdal.GCP(x1, y1, 0, 0, 0),
    gdal.GCP(x2, y2, 0, 1024, 768),
    gdal.GCP(xm1, ym1, 0, 1024, 0),
    gdal.GCP(xm2, ym2, 0, 0, 768),
    gdal.GCP(xc, yc, 0, 512, 384),
    gdal.GCP(xc1, yc1, 0, 512, 0),
    gdal.GCP(xc2, yc2, 0, 0, 384),
    gdal.GCP(xcm1, ycm1, 0, 1024, 384),
    gdal.GCP(xcm2, ycm2, 0, 512, 768)]

ds.SetGCPs(gcps, "EPSG:4326")

geotransform = gdal.GCPsToGeoTransform(gcps)
ds.SetGeoTransform(geotransform)

srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
ds.SetProjection(srs.ExportToWkt())

ds = None

This is the result I am getting enter image description here

This is the result of georefing and mosaicing about 150 images I have total of some 1000+ images

This is what I want enter image description here

Leigh Bettenay
  • 1,401
  • 8
  • 16
Moez
  • 13
  • 4
  • 1
    I know nothing about mars projections, but lat/longs usually have positive/negative numbers, just like in your second picture, so my guess is that your input values perhaps need to be looked at for a negative sign somewhere in there. – nr_aus Mar 14 '23 at 01:14
  • Yeah it does contain negative values, and the georef images does lie on both positive and negative coordinates, – Moez Mar 14 '23 at 17:41
  • I think the format for gdal.GCP, should be 'gdal.GCP(pixel_x, pixel_y, longitude, latitude). have a look at this article which shows how gcp is used there, might be a clue somewhere there https://gis.stackexchange.com/questions/327574/georeference-an-unreferenced-image-using-rasterio-in-python – nr_aus Mar 15 '23 at 01:47

0 Answers0