1

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:

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 gdal
ds = 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 area
data = 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 distance

The 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?

Community
  • 1
TommasoF
  • 143
  • 6

1 Answers1

1

You can simply reproject your shapefile to a projection that has meters as units, like the UTM zone of your part of the world.

Don't use Google/Pseudo/World mercator projection, its units coincide with real meters only at the aequator.

AndreJ
  • 76,698
  • 5
  • 86
  • 162
  • Thank you for your answer. I've reprojected the coordinates to 'WGS84/UTM zone 33N' but after executing gdal_proximity script, the result distance is expressed in degrees... – TommasoF Jul 03 '14 at 07:28
  • Can you add how you did the reprojection? – AndreJ Jul 03 '14 at 07:46
  • I've used the QGIS utility. I've set the Layer CSR and saved it as shapefile. Now the selected CSR is WGS 84 / UTM zone 33N. – TommasoF Jul 03 '14 at 07:53
  • 1
    No, Set Layer CRSis definitively wrong. You corrupt your data with that because it does not change the coordinate values. So change it back, and use Save As... and choose another CRS there. – AndreJ Jul 03 '14 at 08:08
  • Unfortunately every time I tried to follow your instructions, Qgis crashes... There is an other way to do it? Like a gdal script? – TommasoF Jul 03 '14 at 08:42
  • 1
    I was mistaken by the words "shape file" you used. For reprojecting rasters, gdalwarp would be best. QGIS might crash when the raster file is large. Keep in mind that UTM 33N is not vaild worldwide, so you might have to clip the source file. – AndreJ Jul 03 '14 at 09:15