Hello, I am trying to reproject a population distribution raster file of China into 3857 space. I tried but the images are not aligning. Looking at the CRS it has a lot of different standards, I am not sure what to make of it.
The meta of the source file looks something like this:
{'driver': 'GTiff',
'dtype': 'int32',
'nodata': 2147483647.0,
'width': 5049,
'height': 5582,
'count': 1,
'crs': CRS.from_wkt('PROJCS["Albers_Conic_Equal_Area",GEOGCS["GCS_Unknown_datum_based_upon_the_Krassowsky_1940_ellipsoid",DATUM["Not_specified_based_on_Krassowsky_1940_ellipsoid",SPHEROID["Krasovsky_1940",6378245,298.3,AUTHORITY["EPSG","7024"]],AUTHORITY["EPSG","6024"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["standard_parallel_1",25],PARAMETER["standard_parallel_2",47],PARAMETER["latitude_of_center",0],PARAMETER["longitude_of_center",105],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]]]'),
'transform': Affine(1000.0, 0.0, -2638055.0614613006,
0.0, -1000.0, 5922165.484718928)}
PopTIF.crs.data returns:
{'proj': 'aea',
'lat_1': 25,
'lat_2': 47,
'lat_0': 0,
'lon_0': 105,
'x_0': 0,
'y_0': 0,
'ellps': 'krass',
'units': 'm',
'no_defs': True}
after running the reprojection algorithm:
{'driver': 'GTiff',
'dtype': 'int32',
'nodata': 2147483647.0,
'width': 6708,
'height': 6102,
'count': 1,
'crs': CRS.from_dict(init='epsg:3857'),
'transform': Affine(1184.7466003556895, 0.0, 7551104.353792087,
0.0, -1184.7466003556895, 7299154.069454485)}
I am then using PopTIF.dataset_mask() to mask the area of a Heightfield GEOTIFF (I requested the raster area of this file to be the same as the lat/lon of the PopTIF data:
'lat_1': 25,'lat_2': 47, 'lat_0': 0, 'lon_0': 105
{'driver': 'GTiff',
'dtype': 'float32',
'nodata': -32767.0,
'width': 2816,
'height': 2288,
'count': 1,
'crs': CRS.from_dict(init='epsg:4326'),
'transform': Affine(0.02197265625, 0.0, 73.4765625,
0.0, -0.02197265625, 53.7890625)}
EPSG:4326 to EPSG:3857
{'driver': 'GTiff',
'dtype': 'float32',
'nodata': -32767.0,
'width': 2594,
'height': 2537,
'count': 1,
'crs': CRS.from_dict(init='epsg:3857'),
'transform': Affine(2655.775407243914, 0.0, 8179373.5227401415,
0.0, -2655.775407243914, 7130308.041557059)}
So, After moving both GEOTIFFs into 3857 space. I then needed to upscale the Heightfield so the array matched the the size of the Population data:
import skimage.transform as st
pathToRaster = r'./DEMTIF/cnl10_3857.tif'
newsize = reprojectedmask.shape
with rio.open(pathToRaster) as src:
arr = src.read()
arr[arr ==-9999] = np.nan
arr = arr[0,:,:]
new_shape = newsize
newarr= st.resize(arr, new_shape, mode='constant')
Resulting array can be masked, and looks something like this, Notice how there is no height data in areas such as Taiwan, There has been a mismatch with the projections, However I cant tell where/how. This may be due to the fact I specified the lat/lon in AEA(plus other CRS)? Would appreciate any guidance on how to find the correct area to sample and how to transform it
