Last active
October 3, 2021 19:17
-
-
Save scidam/9727937cdd713182f1168be2a42eaeb4 to your computer and use it in GitHub Desktop.
Getting area from tiffs file (SDM model)
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| # -*- coding: utf-8 -*- | |
| """ | |
| Created on Sun Jun 7 07:07:17 2020 | |
| @author: dmitry | |
| """ | |
| from osgeo import gdal | |
| import numpy as np | |
| def get_area(filename='',maskfile='', third=None, threshold=0.8): | |
| R = 6631000.0 | |
| ds = gdal.Open(filename, gdal.GA_ReadOnly) | |
| rb = ds.GetRasterBand(1) | |
| img_array = rb.ReadAsArray() | |
| msk = gdal.Open(maskfile, gdal.GA_ReadOnly) | |
| rbm = msk.GetRasterBand(1) | |
| img_mask = rbm.ReadAsArray() | |
| gt=ds.GetGeoTransform() | |
| dlat = gt[-1] | |
| dlon = gt[1] | |
| m, n = img_array.shape | |
| lats = np.linspace(gt[-3], gt[-3] + dlat * n, n) | |
| lons = np.linspace(gt[0], gt[0] + dlon * m, m) | |
| LON, LAT = np.meshgrid(lons, lats) | |
| if third: | |
| msk = gdal.Open(third, gdal.GA_ReadOnly) | |
| rbm = msk.GetRasterBand(1) | |
| third = rbm.ReadAsArray() | |
| with_mask = (R*np.cos(LAT[(img_array > threshold)*(img_mask > threshold)*(third > threshold)]/180*np.pi)*R*np.abs(dlat/180*dlon/180*np.pi*np.pi)).sum()/10**6 | |
| else: | |
| with_mask = (R*np.cos(LAT[(img_array > threshold)*(img_mask > threshold)]/180*np.pi)*R*np.abs(dlat/180*dlon/180*np.pi*np.pi)).sum()/10**6 | |
| without_mask =(R*np.cos(LAT[img_array > threshold]/180*np.pi)*R*np.abs(dlat/180*dlon/180*np.pi*np.pi)).sum()/10**6 | |
| return with_mask, without_mask | |
| filenames = ['pinus__70cc26_RF_100_0'] | |
| computing_plan = [ | |
| # ('picea__70mr85_RF_100_0.tiff', 'picea_0.tiff', 0.434), | |
| # ('picea__70mr26_RF_100_0.tiff', 'picea_0.tiff', 0.434), | |
| # ('abies nephrolepis__70mr85_RF_100_0.tiff', 'abies nephrolepis_0.tiff', 0.404), | |
| # ('abies nephrolepis__70mr26_RF_100_0.tiff', 'abies nephrolepis_0.tiff', 0.404), | |
| # ('abies sachalinensis__70mr85_RF_100_0.tiff', 'abies sachalinensis_0.tiff', 0.404), | |
| # ('abies sachalinensis__70mr26_RF_100_0.tiff', 'abies sachalinensis_0.tiff', 0.404), | |
| # ('abies__70mr85_RF_100_0.tiff', 'abies_0.tiff', 0.414), | |
| # ('abies__70mr26_RF_100_0.tiff', 'abies_0.tiff', 0.414), | |
| # ('picea__mrlgm_RF_100_0.tiff', 'picea_0.tiff', 0.434), | |
| # ('abies nephrolepis__mrlgm_RF_100_0.tiff', 'abies nephrolepis_0.tiff', 0.404), | |
| # ('abies sachalinensis__mrlgm_RF_100_0.tiff', 'abies sachalinensis_0.tiff', 0.404), | |
| # ('abies__mrlgm_RF_100_0.tiff', 'abies_0.tiff', 0.414) | |
| # ('picea_0.tiff', 'picea_0.tiff', 0.434), | |
| # ('abies nephrolepis_0.tiff', 'abies nephrolepis_0.tiff', 0.404), | |
| # ('abies sachalinensis_0.tiff', 'abies sachalinensis_0.tiff', 0.404), | |
| # ('abies_0.tiff', 'abies_0.tiff', 0.414) | |
| ('pinus_0.tiff', 'pinus__mrlgm_RF_100_0.tiff', 'pinus__70mr85_RF_100_0.tiff', 0.2), | |
| ('pinus_0.tiff', 'pinus__mrlgm_RF_100_0.tiff', 'pinus__70mr26_RF_100_0.tiff', 0.2) | |
| ] | |
| for name, mask, tr, th in computing_plan: | |
| area_masked, area = get_area(filename=name, maskfile=mask, third=tr, threshold=th) | |
| print(f"S w/ mask = {area_masked} sq. km; S w/o mask={area} sq. km; filename={name}, threshold={th}. maskfile = {mask},{tr}.") |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment