fork(7) download
copy 
  1. import numpy as np
  2. import gdal
  3. from osgeo import gdal
  4. from osgeo import osr
  5. from osgeo import ogr
  6. from osgeo.gdalconst import *
  7. gdal.AllRegister() # register all drivers
  8. gdal.UseExceptions()
  9.  
  10. '''http://m...content-available-to-author-only...l.com/blog/archive/2012/5/2/understanding-raster-basic-gis-concepts-and-the-python-gdal-library/'''
  11.  
  12. #############
  13. # Functions #
  14. #############
  15.  
  16. def transform_utm_to_wgs84(easting, northing, zone):
  17.     utm_coordinate_system = osr.SpatialReference()
  18.  
  19.     # Set geographic coordinate system to handle lat/lon
  20.     utm_coordinate_system.SetWellKnownGeogCS("WGS84") 
  21.     is_northern = northing > 0    
  22.     utm_coordinate_system.SetUTM(zone, is_northern)
  23.  
  24.     # Clone ONLY the geographic coordinate system 
  25.     wgs84_coordinate_system = utm_coordinate_system.CloneGeogCS() 
  26.  
  27.     # create transform component
  28.     utm_to_wgs84_geo_transform = osr.CoordinateTransformation(utm_coordinate_system, wgs84_coordinate_system) # (, )
  29.  
  30.     # returns lon, lat, altitude
  31.     return utm_to_wgs84_geo_transform.TransformPoint(easting, northing, 0) 
  32.  
  33. class WGS84Transform(object):
  34.     # TODO decide whether to have init and/or call functions
  35.     #def transform_wgs84_to_utm(lon, lat)
  36.  
  37. #    def __init__(self,lon,lat):
  38. #    self.lon = lon
  39. #	self.lat = lat
  40.  
  41.     def get_utm_zone(self,longitude):
  42. 	return (int(1+(longitude+180.0)/6.0))
  43.  
  44.     def is_lat_northern(self,latitude):
  45. 	"""
  46. 	Determines if given latitude is a northern for UTM
  47. 	"""
  48. 	if (latitude < 0.0):
  49. 	    return 0
  50. 	else:
  51. 	    return 1
  52.  
  53.     def wgs84_to_utm(self,lon,lat):            
  54. 	    utm_coordinate_system = osr.SpatialReference()
  55. 	    # Set geographic coordinate system to handle lat/lon  
  56. 	    utm_coordinate_system.SetWellKnownGeogCS("WGS84") 
  57. 	    utm_coordinate_system.SetUTM(self.get_utm_zone(lon), self.is_lat_northern(lat))
  58.  
  59. 	    # Clone ONLY the geographic coordinate system  
  60. 	    wgs84_coordinate_system = utm_coordinate_system.CloneGeogCS() 
  61.  
  62. 	    # create transform component
  63. 	    wgs84_to_utm_geo_transform = osr.CoordinateTransformation(wgs84_coordinate_system, utm_coordinate_system) # (, )
  64. 	    # returns easting, northing, altitude
  65. 	    return wgs84_to_utm_geo_transform.TransformPoint(lon, lat, 0) 
  66.  
  67. def get_iterable_extent(*args):
  68.     '''Returns list of minimum and maximum from lists/array input'''
  69.     iterable_extent = list()
  70.     for iter_object in args:
  71. 	iterable_extent.append(min(iter_object))
  72. 	iterable_extent.append(max(iter_object))
  73.     return iterable_extent 
  74.  
  75.  
  76. def get_raster_size(min_x, min_y, max_x, max_y, cell_width, cell_height):
  77.     """Determine the number of rows/columns given the bounds of the point 
  78.     data and the desired cell size"""
  79.  
  80.     print 'raster min_x:',min_x 
  81.     print 'raster max_x:',max_x 
  82.     print 'raster min_y:',min_y 
  83.     print 'raster max_y:',max_y 
  84.     cols = int((max_x - min_x) / cell_width)
  85.     rows = int((max_y - min_y) / abs(cell_height))
  86.     return cols, rows
  87.  
  88.  
  89. def lonlat_to_pixel(lon, lat, inverse_geo_transform):
  90.     """Translates the given lon, lat to the grid pixel coordinates
  91.     in data array (zero start)"""
  92.  
  93.     wgs84_object = WGS84Transform()
  94.     # transform to utm
  95.     utm_x, utm_y, utm_z = wgs84_object.wgs84_to_utm(lon, lat)
  96.     print 'utm_x:',utm_x
  97.     print 'utm_y:',utm_y
  98.     print 'utm_z:',utm_z
  99.  
  100.     # apply inverse geo tranform
  101.     pixel_x, pixel_y = gdal.ApplyGeoTransform(inverse_geo_transform, utm_x, utm_y)
  102.     print 'pixel_x:',pixel_x
  103.     print 'pixel_y:',pixel_y
  104.     pixel_x = int(pixel_x) - 1 # adjust to 0 start for array
  105.     pixel_y = int(pixel_y) - 1 # adjust to 0 start for array
  106.  
  107.  
  108.     return pixel_x, abs(pixel_y) # y pixel is likly a negative value given geo_transform
  109.  
  110.  
  111. def create_raster(lons,lats,values,filename="test.tiff",output_format="GTiff"):      
  112.     """lon/lat values in WGS84"""          
  113.  
  114.     # create empty raster
  115.     driver = gdal.GetDriverByName(output_format)
  116.     number_of_bands = 1
  117.     band_type = gdal.GDT_Float32
  118.  
  119.     x_rotation = 0
  120.     y_rotation = 0
  121.     cell_width_meters = 50.0
  122.     cell_height_meters = 50.0
  123.  
  124.     # retrieve bounds for point data           
  125.     min_lon, max_lon, min_lat, max_lat = get_iterable_extent(lons,lats) 
  126.     print 'min_lon:',min_lon 
  127.     print 'max_lon:',max_lon 
  128.     print 'min_lat:',min_lat 
  129.     print 'max_lat:',max_lat 
  130.  
  131.     # Set geographic coordinate system to handle lat/lon        
  132.     srs = osr.SpatialReference()
  133.     srs.SetWellKnownGeogCS("WGS84")
  134.  
  135.     # Set projected coordinate system  to handle meters         
  136.     wgs84_obj = WGS84Transform()
  137.     srs.SetUTM(wgs84_obj.get_utm_zone(min_lon), wgs84_obj.is_lat_northern(max_lat)) 
  138.  
  139.     # create transforms for point conversion
  140.     wgs84_coordinate_system = srs.CloneGeogCS() # clone only the geographic coordinate system   
  141.     wgs84_to_utm_transform = osr.CoordinateTransformation(wgs84_coordinate_system, srs)
  142.     utm_to_wgs84_transform = osr.CoordinateTransformation(srs, wgs84_coordinate_system)
  143.  
  144.     # convert to UTM
  145.     top_left_x, top_left_y, z = wgs84_obj.wgs84_to_utm(min_lon, max_lat)    
  146.     lower_right_x, lower_right_y, z = wgs84_obj.wgs84_to_utm(max_lon, min_lat)
  147.     print 'top_left_x:',top_left_x
  148.     print 'top_left_y:',top_left_y
  149.     print 'lower_right_x:',lower_right_x
  150.     print 'lower_right_y:',lower_right_y
  151.  
  152.     '''get_raster_size(min_x, min_y, max_x, max_y, cell_width, cell_height)'''
  153.     cols, rows = get_raster_size(top_left_x,lower_right_y,
  154.                                  lower_right_x,top_left_y,
  155.                                  cell_width_meters,cell_height_meters)
  156.     print cols, rows 
  157.     dataset = driver.Create(filename, cols, rows, number_of_bands, band_type)
  158.  
  159.     # GeoTransform parameters
  160.     # --> need to know the area that will be covered to define the geo tranform
  161.     # top left x, w-e pixel resolution, rotation, top left y, rotation, n-s pixel resolution
  162.     geo_transform = [ top_left_x, cell_width_meters, x_rotation, top_left_y, y_rotation, -cell_height_meters ] # cell height must be negative (-) to apply image space to map
  163.     dataset.SetGeoTransform(geo_transform)
  164.     dataset.SetProjection(srs.ExportToWkt())
  165.     inverse_geo_transform = gdal.InvGeoTransform(geo_transform)[1] # for mapping lat/lon to pixel
  166.  
  167.     # get the empty raster data array
  168.     band = dataset.GetRasterBand(1) # 1 == band index value
  169.     data = band.ReadAsArray(0, 0, cols, rows).astype(np.cfloat)
  170.  
  171.     # TODO check that lat/lon don't need to get meshed
  172.     # populate array values for output
  173.     for lon, lat, value in zip(lons,lats,values):
  174.         # apply value to array
  175. 	#TODO figure out why pixel_x is negative
  176.         pixel_x, pixel_y = lonlat_to_pixel(lon, lat, inverse_geo_transform)
  177. 	print pixel_x, pixel_y, data.shape
  178.         data[pixel_x][pixel_y] = value
  179.  
  180.     # write the updated data array to file
  181.     band.WriteArray(data, 0, 0)
  182.     band.SetNoDataValue(NULL_VALUE)
  183.     band.FlushCache()
  184.  
  185.     # set dataset to None to "close" file
  186.     dataset = None
  187.     # TODO have function return pixel values for interpolation
  188.     return pixel_x, pixel_y
  189.  
  190. #################
  191. # Main Function #
  192. #################
  193.  
  194. if __name__ == '__main__':
  195. 	# example coordinates
  196. 	lat = [45.3,56.2,23.4]
  197. 	lon = [134.6,128.7,111.9]
  198. 	val = [3,6,2]
  199.  
  200. 	create_raster(lon,lat,val)
Runtime error #stdin #stdout 0.08s 8832KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
Standard output is empty
https://ideone.com/CusKmZ
language:
Python (cpython 2.7.16)
created:
11 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!