使用rasterio在python中执行栅格计算

iovurdzv  于 2023-08-02  发布在  Python
关注(0)|答案(1)|浏览(159)

我尝试将两个光栅(raster1和raster2)相乘,然后使用相乘的结果写入光栅。
我可以使用python中的rasterio包读取输入光栅,但我很难将输出光栅写入文件。
我哪里做错了?

  1. rasterio version 1.3.8
  2. Python 3.10.8
  3. The meta data for the input rasters are:
  4. raster1:
  5. Name JRC_forest_yield_upsampled_to_1km
  6. Path C:\Users\Input_datasets\Forest_Biomass_Map\JRC_forest_yield_upsampled_to_1km.tif
  7. Extent 900000.0000000000000000,900000.0000000000000000 : 7400000.0000000000000000,5416000.0000000000000000
  8. Width 6500
  9. Height 4516
  10. Data type Float64 - Sixty four bit floating point
  11. GDAL Driver Description GTiff
  12. GDAL Driver Metadata GeoTIFF
  13. Compression LZW
  14. CRS EPSG:3035 - ETRS89-extended / LAEA Europe
  15. raster2:
  16. Name eu_nuts0_for_rents_prices_rast
  17. Path C:\Users\Throughput_datasets\eu_nuts0_for_rents_prices_rast.tif
  18. Extent 900000.0000000000000000,900000.0000000000000000 : 7400000.0000000000000000,5416000.0000000000000000
  19. Width 6500
  20. Height 4516
  21. Data type Float32 - Thirty two bit floating point
  22. GDAL Driver Description GTiff
  23. GDAL Driver Metadata GeoTIFF
  24. Compression (blank in QGIS)
  25. CRS EPSG:3035 - ETRS89-extended / LAEA Europe

字符串
错误代码:

  1. (.venv) C:\Users\spencerd\Documents> raster_calculator.py
  2. c:\Users\spencerd\Documents\.venv\lib\site-packages\rasterio\__init__.py:329: NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix will be returned.
  3. dataset = writer(


脚本:

  1. import rasterio
  2. from rasterio import features
  3. from rasterio.enums import MergeAlg
  4. from rasterio.plot import show
  5. from numpy import int16
  6. # Paths to the input raster files
  7. raster1_path = 'C:/Users/Spatial datasets/EULandClass/EU_landSystem.tif'
  8. raster_woody_path = 'C:/Users/Spatial datasets/European countries/Input_datasets/Forest_Biomass_Map/JRC_forest_yield_upsampled_to_1km.tif'
  9. raster2_path = 'C:/Users/Throughput_datasets/eu_nuts0_for_rents_prices_rast.tif'
  10. # Open the input raster files and plot
  11. rast1 = rasterio.open(raster_woody_path).read()
  12. rast2 = rasterio.open(raster2_path).read()
  13. result = rast1 * rast2
  14. result
  15. # define output file path and parameters
  16. output_path = "C:/Users/Throughput_datasets/eu_forestry_monetary_yields_unmasked.tif"
  17. # Write the result to a new raster file
  18. with rasterio.open(output_path,
  19. 'w',
  20. width=6500,
  21. height=4516,
  22. count=1,
  23. crs='epsg:3050',
  24. dtype='float32',
  25. nodata=-0,
  26. driver='GTiff',
  27. transform= None) as dst:
  28. dst.write(result)


从这段代码中,我想将结果光栅“result”保存到文件路径“output_path”,但我只是得到一个错误,我需要做什么不同的事情?

n3ipq98p

n3ipq98p1#

最简单的方法是从一个用作源的栅格中获取所需的参数(包括宽度,高度,crs等)。
按以下方式编辑最后一段代码:

  1. # Write the result to a new raster file
  2. kwargs = rast1.meta
  3. kwargs.update(
  4. dtype=rasterio.float32,
  5. count=1)
  6. with rasterio.open(output_path, 'w', **kwargs) as dst:
  7. dst.write_band(1, result.astype(rasterio.float32))

字符串

相关问题