# Get pixel(Row Column) from Lat Long for Moon

How to obtain pixel(Row Column) Value from Latitude and Longitude given, for Moon?

Actually, I am Having a GeoTiff file which contains craters, and I need to crop the specific Crater whose Center is given to me in the form of Latitude and longitude.

I know that I have to use gdal to get this done, but I am having wrong results from my code.

Code is in Python Code:

``````from osgeo import osr, gdal
import sys

# get the existing coordinate system
# create the new coordinate system
ds = gdal.Open("C:\Users\ZARL\Desktop\PRL_Project\small.tif")
old_cs= osr.SpatialReference()
old_cs.ImportFromWkt(ds.GetProjectionRef())
wgs84_wkt = """
GEOGCS["GCS_MOON",
DATUM["D_MOON",
PRIMEM["Reference_Meridian",0],
UNIT["degree",0.0174532925199433]],
PROJECTION["Equirectangular"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",22.46049118041992],
PARAMETER["standard_parallel_1",0],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]]]"""
new_cs = osr.SpatialReference()
new_cs .ImportFromWkt(wgs84_wkt)

# create a transform object to convert between coordinate systems
transform = osr.CoordinateTransformation(old_cs,new_cs)

#get the point to transform, pixel (0,0) in this case
width = ds.RasterXSize
height = ds.RasterYSize
gt = ds.GetGeoTransform()
minx = gt
miny = gt + width*gt + height*gt

#get the coordinates in lat long
#latlong = transform.TransformPoint(minx,miny)

transform1 = ds.GetGeoTransform()

xOrigin = transform1
yOrigin = transform1
pixelWidth = transform1
pixelHeight = -transform1

latlong = transform.TransformPoint(xOrigin,yOrigin)
#data = band.ReadAsArray(0, 0, width, height)

points_list = [ (22.4655,-71.7136)] #list of Lat long

for point in points_list:
col = int((point - latlong) / pixelWidth)
row = int((latlong-point ) / pixelHeight)

print row,col
``````

I am getting wrong values.

The Tiff file is:https://www.dropbox.com/s/twov5jpe40hcf4q/small.tif?dl=0

when I run gdalinfo command on the picture I get the following:

``````    Driver: GTiff/GeoTIFF
Files: small.tif
small.tif.aux.xml
Size is 40, 33
Coordinate System is:
PROJCS["unnamed",
GEOGCS["GCS_MOON",
DATUM["D_MOON",
PRIMEM["Reference_Meridian",0],
UNIT["degree",0.0174532925199433]],
PROJECTION["Equirectangular"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",22.46049118041992],
PARAMETER["standard_parallel_1",0],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]]]
Origin = (135.724395353784870,-2174581.951404873300000)
Pixel Size = (0.812720930262190,-0.812720930262190)
AREA_OR_POINT=Area
TIFFTAG_XRESOLUTION=1
TIFFTAG_YRESOLUTION=1
INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (     135.724,-2174581.951) ( 22d27'53.88"E, 71d42'47.22"S)
Lower Left  (     135.724,-2174608.771) ( 22d27'53.88"E, 71d42'50.40"S)
Upper Right (     168.233,-2174581.951) ( 22d27'57.74"E, 71d42'47.22"S)
Lower Right (     168.233,-2174608.771) ( 22d27'57.74"E, 71d42'50.40"S)
Center      (     151.979,-2174595.361) ( 22d27'55.81"E, 71d42'48.81"S)
Band 1 Block=40x1 Type=Byte, ColorInterp=Gray
Min=1.000 Max=254.000
Minimum=1.000, Maximum=254.000, Mean=122.977, StdDev=48.776
STATISTICS_MAXIMUM=254
STATISTICS_MEAN=122.97651515152
STATISTICS_MINIMUM=1
STATISTICS_STDDEV=48.776222721807
``````

Any Other Approach is also fine.

• Is the WKT working at all? It should start out "PROJCS["name",GEOGCS{]..." – mkennedy Feb 7 '18 at 18:26
• That "Tiff" is a PNG and without the original Tiff we can't get the metadata that it might contain. For example, we need to know the corner coordinates of the image, which are normally obtained from GeoTIFF metadata. – Spacedman Feb 7 '18 at 18:53
• I have updated the question – abhinav antani Feb 9 '18 at 5:29

I have done the solution and it is working for almost every cases. I am putting the solution over here. It is done in python.

Code:

``````from osgeo import osr,gdal

infile = "C:\Users\ZARL\Desktop\PRL_Project\gale463m.tif" #File Path goes here
#Lat Long value
long = 137.855486
lat =  -5.378912

indataset = gdal.Open( infile)

srs = osr.SpatialReference()
srs.ImportFromWkt(indataset.GetProjection())

srsLatLong = srs.CloneGeogCS()
ct = osr.CoordinateTransformation(srsLatLong, srs)
(X, Y, height) = ct.TransformPoint(long, lat)

# Report results
print('longitude: %f\t\tlatitude: %f' % (long, lat))
print('X: %f\t\tY: %f' % (X, Y))
#VALUE OF COORDINATE IN METERS
print X ,Y

driver = gdal.GetDriverByName('GTiff')
band = indataset.GetRasterBand(1)

cols = indataset.RasterXSize
rows = indataset.RasterYSize

transform = indataset.GetGeoTransform()

xOrigin = transform
yOrigin = transform
pixelWidth = transform
pixelHeight = -transform
data = band.ReadAsArray(0, 0, cols, rows)

points_list = [(X,Y)] #list of X,Y coordinates

for point in points_list:
col = int((point - xOrigin) / pixelWidth )
row = int((yOrigin - point ) / pixelHeight)
#ROW AND COLUMN VALUE
print row,col
#Data AT THAT ROW COLUMN
value = data[row][col]

print value
``````

This will give the data value at that pixel whose Lat Long is provided.