I followed the answer given in this topic for searching the distance to the nearest coastline given a lat/lon coordinate. The distance output obtained is in degree (I think) but I need the distance expressed in meters.
By citing the answer given:
from osgeo import gdalds = gdal.Open("/home/user/distance.tif")gt = ds.GetGeoTransform()topLeftX = gt[0]dimX = gt[1]topLeftY = gt[3]dimY = gt[5]band = ds.GetRasterBand(1)#Numpy array of the distances in the whole areadata = band.ReadAsArray()#Indexing by longitude/latitude position def getDistance(coord): x = int((coord[0] - topLeftX)/dimX) y = int((coord[1] - topLeftY)/dimY) distance = data[y,x] return distanceThe coordinate system is WGS84. I think I have to rasterize the shape file in a Projected Coordinate System but then I don't know what to do. Any suggestions?
أكثر...
By citing the answer given:
There will be a units problem if you create the distance raster in geographic coordinates. In order to get kilimetres you will need to rasterize the sea in a suitable projected coordinate system and then project the distance raster to geograpic coordinates to do the lookup. This will give you the fastest response time as you're not projecting the query coordinate at the time of lookup.
At the moment my python code is the following:
from osgeo import gdalds = gdal.Open("/home/user/distance.tif")gt = ds.GetGeoTransform()topLeftX = gt[0]dimX = gt[1]topLeftY = gt[3]dimY = gt[5]band = ds.GetRasterBand(1)#Numpy array of the distances in the whole areadata = band.ReadAsArray()#Indexing by longitude/latitude position def getDistance(coord): x = int((coord[0] - topLeftX)/dimX) y = int((coord[1] - topLeftY)/dimY) distance = data[y,x] return distanceThe coordinate system is WGS84. I think I have to rasterize the shape file in a Projected Coordinate System but then I don't know what to do. Any suggestions?
أكثر...