I want to create an "empty" raster with predefined coordinates from the upper left corner and a given Pixelsize. My first step is to calculate the rows, and columns. Therefore the distance between the upper left and upper right corner and the upper left and lower left corner will be calculated by the given code:
upper_left_latlon = [52.00, 9.50]lower_right_lat_lon = [50.00, 13.00]from math import sin, cos, sqrt, atan2, radiansdef coord_dist(lat1, long1, lat2, long2): # source: NSSDC - Earth Fact Sheet r_earth = 6378.1 dist_long = radians(long2) - radians(long1) dist_lat = radians(lat2) - radians(lat1) a = sin(dist_lat / 2)**2 + cos(lat1) * cos(lat2) * sin(dist_long / 2)**2 c = 2 * atan2(sqrt(a),sqrt(1 - a)) return((r_earth*c)*1000)col = roundup(dist_latlong.coord_dist(upper_left_latlon[0],upper_left_latlon[1],upper_left_latlon[0],lower_right_lat_lon[1])/size)row = roundup(dist_latlong.coord_dist(upper_left_latlon[0],upper_left_latlon[1],lower_right_lat_lon[0],upper_left_latlon[1])/size)After this calculation i want to create a raster with gdal. The result has the correct pixelsize but something went wrong with the projection. I think there is a lapse of thought in my code by transforming the raster.
driver = gdal.GetDriverByName('GTiff')outRaster = driver.Create('test_raster.tif',col,row,1,gdal.GDT_Int32)outRaster.SetGeoTransform((upper_left_latlon[0],size,0,upper_left_latlon[1],0,-size))outband = outRaster.GetRasterBand(1)outRasterSRS = osr.SpatialReference()outRasterSRS.ImportFromEPSG(4326)outRaster.SetProjection(outRasterSRS.ExportToWkt())outband.FlushCache()
أكثر...
upper_left_latlon = [52.00, 9.50]lower_right_lat_lon = [50.00, 13.00]from math import sin, cos, sqrt, atan2, radiansdef coord_dist(lat1, long1, lat2, long2): # source: NSSDC - Earth Fact Sheet r_earth = 6378.1 dist_long = radians(long2) - radians(long1) dist_lat = radians(lat2) - radians(lat1) a = sin(dist_lat / 2)**2 + cos(lat1) * cos(lat2) * sin(dist_long / 2)**2 c = 2 * atan2(sqrt(a),sqrt(1 - a)) return((r_earth*c)*1000)col = roundup(dist_latlong.coord_dist(upper_left_latlon[0],upper_left_latlon[1],upper_left_latlon[0],lower_right_lat_lon[1])/size)row = roundup(dist_latlong.coord_dist(upper_left_latlon[0],upper_left_latlon[1],lower_right_lat_lon[0],upper_left_latlon[1])/size)After this calculation i want to create a raster with gdal. The result has the correct pixelsize but something went wrong with the projection. I think there is a lapse of thought in my code by transforming the raster.
driver = gdal.GetDriverByName('GTiff')outRaster = driver.Create('test_raster.tif',col,row,1,gdal.GDT_Int32)outRaster.SetGeoTransform((upper_left_latlon[0],size,0,upper_left_latlon[1],0,-size))outband = outRaster.GetRasterBand(1)outRasterSRS = osr.SpatialReference()outRasterSRS.ImportFromEPSG(4326)outRaster.SetProjection(outRasterSRS.ExportToWkt())outband.FlushCache()
أكثر...