SHARE
TWEET

Untitled

a guest Feb 22nd, 2019 68 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. ### Feng Yin
  2. ### Department of Geography, UCL
  3. ### ucfafyi@ucl.ac.uk
  4.  
  5. import osr          
  6. import ogr          
  7. import gdal        
  8. import json        
  9. import numpy  as np
  10.                    
  11. def get_boundary(raster, to_wgs84 = True):
  12.     '''            
  13.     Take a raster and find boundary of the raster
  14.     if to_wgs84 = True, transform them into wgs84 coordinates
  15.     with at least meters resoluton.
  16.     Return geojson(s)                                                                                                                                                                
  17.     '''            
  18.     g    = gdal.Open(raster)
  19.     geom = g.GetGeoTransform()
  20.     x_coords = geom[0] + np.arange(g.RasterXSize) * geom[1]
  21.     y_coords = geom[3] + np.arange(g.RasterYSize) * geom[5]
  22.     north_boundary = np.dstack([x_coords,       np.ones(g.RasterXSize) * y_coords[0]])
  23.     south_boundary = np.dstack([x_coords[::-1], np.ones(g.RasterXSize) * y_coords[-1]])
  24.     east_boundary  = np.dstack([x_coords[-1] * np.ones(g.RasterYSize), y_coords])
  25.     west_boundary  = np.dstack([x_coords[0]  * np.ones(g.RasterYSize), y_coords[::-1]])
  26.     boundary = np.hstack([north_boundary, east_boundary, south_boundary, west_boundary]).tolist()
  27.     polygon  = {'type': 'Polygon', 'coordinates': boundary}
  28.     geometry = ogr.CreateGeometryFromJson(json.dumps(polygon))
  29.     sgeometry = geometry.SimplifyPreserveTopology(abs(geom[1])/1000)
  30.     spolygon = json.loads(sgeometry.ExportToJson())
  31.     sourceprj = osr.SpatialReference()
  32.     sourceprj.ImportFromWkt(g.GetProjection())
  33.     try:            
  34.         urn = sourceprj.ExportToXML().split('srsID')[1].split('=')[1].split('</gml:name>')[0].replace('"', '').replace('>', '')
  35.         geojson = {
  36.           "type": "FeatureCollection",
  37.           "name": "AOI",
  38.           "crs": { "type": "name", "properties": { "name": urn} },
  39.           "features": [
  40.           { "type": "Feature", "properties": {}, "geometry": spolygon}
  41.                       ]
  42.          }          
  43.     except:        
  44.         proj4 = sourceprj.ExportToProj4()
  45.         geojson = {
  46.           "type": "FeatureCollection",
  47.           "name": "AOI",
  48.           "crs": { "type": "name", "properties": { "name": proj4, "type": "proj4"} },
  49.           "features": [
  50.           { "type": "Feature", "properties": {}, "geometry": spolygon}
  51.                       ]
  52.          }          
  53.                    
  54.     geojson = json.dumps(geojson)
  55.                    
  56.     if to_wgs84:    
  57.         targetprj = osr.SpatialReference()
  58.         targetprj.ImportFromEPSG(4326)
  59.         transform = osr.CoordinateTransformation(sourceprj, targetprj)
  60.         geometry.Transform(transform)
  61.                    
  62.         stgeometry = geometry.SimplifyPreserveTopology(abs(geom[1]) / 10. * 0.000001)
  63.         stpolygon = json.loads(stgeometry.ExportToJson())
  64.                    
  65.         tgeojson = {
  66.           "type": "FeatureCollection",
  67.           "name": "AOI",
  68.           "crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:EPSG::4326"} },
  69.           "features": [
  70.           { "type": "Feature", "properties": {}, "geometry": stpolygon}
  71.                       ]
  72.          }          
  73.         tgeojson = json.dumps(tgeojson)
  74.                    
  75.         return tgeojson, geojson
  76.     else:          
  77.         return geojson
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top