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
#SPHEROID["MOON_localRadius",1737400,0]
# 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",
SPHEROID["MOON_localRadius",1737400,0]],
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[0]
miny = gt[3] + width*gt[4] + height*gt[5]
#get the coordinates in lat long
#latlong = transform.TransformPoint(minx,miny)
transform1 = ds.GetGeoTransform()
xOrigin = transform1[0]
yOrigin = transform1[3]
pixelWidth = transform1[1]
pixelHeight = -transform1[5]
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[0] - latlong[1]) / pixelWidth)
row = int((latlong[0]-point[1] ) / 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",
SPHEROID["MOON_localRadius",1737400,0]],
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)
Metadata:
AREA_OR_POINT=Area
TIFFTAG_XRESOLUTION=1
TIFFTAG_YRESOLUTION=1
Image Structure Metadata:
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
Metadata:
STATISTICS_MAXIMUM=254
STATISTICS_MEAN=122.97651515152
STATISTICS_MINIMUM=1
STATISTICS_STDDEV=48.776222721807
Any Other Approach is also fine.
Best Answer
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:
This will give the data value at that pixel whose Lat Long is provided.