How to spatially distribute the values of a grid into another grid with different projection and resolution?
Added by Natalia L about 5 years ago
I have an emission database (emission.nc) in WGS84 projection and 0.1x0.1degree. I want to reproject the emission database to an LCC projection with 2km resolution (lcc_grid.txt). So, I need to spatially distribute the values of each grid cell of emission.nc into the grid cells of the lcc grid taking into account the area of the emissions.nc cells that overlaps the lcc grid.
When I use remapycon/remapnn/remapcon it seems that the distribution is not correct as probably if one grid cell of the database has the value of 100kg and overlaps 5cells of the lcc grid then all the cells of lcc takes the value of ~100.
Also, with remapcon I get some negative values. In case of using a land mask, remapcon gives zero values for all the grid.
I work as following
cat > $Database_outdir/grid_reproj1.txt <<EOF
gridtype = projection
gridsize = $GRDSIZE1
xsize = $dx1
ysize = $dy1
xname = x
xlongname = "x coordinate of projection"
xunits = "m"
yname = y
ylongname = "y coordinate of projection"
yunits = "m"
xfirst = $ORIGX
xinc = $RESX1
yfirst = $ORIGY
yinc = $RESY1
grid_mapping = lambert_conformal_conic
grid_mapping_name = lambert_conformal_conic
longitude_of_central_meridian = $central_lon
false_easting = 0.
false_northing = 0.
latitude_of_projection_origin = $central_lat
standard_parallel = $standard_parallel_1 $standard_parallel_2
longitude_of_prime_meridian = 0.
earth_radius=6370000.
inverse_flattening = 0.
longitudeOfFirstGridPointInDegrees = $orig_lon
latitudeOfFirstGridPointInDegrees =$orig_lat
earth_radius = 6370000.
inverse_flattening = 0.
spatial_ref = "PROJCS["unnamed",GEOGCS["unnamed ellipse",DATUM["unknown",SPHEROID["unnamed",6370000, 6370000]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["standard_parallel_1",$standard_parallel_1],PARAMETER["standard_parallel_2",$standard_parallel_2],PARAMETER["latitude_of_origin",$central_lat],PARAMETER["central_meridian",$central_lon],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]]]"
GeoTransform = "747678 $RESX1 0 -1438342 0 -$RESY1 "
EOF
cdo remapcon,$Database_outdir/"lcc_grid.txt" $Database_outdir/emissions.nc $Database_outdir/"tmp_nc"
cdo setgridtype,curvilinear -setgrid,$Database_outdir/"lcc_grid.txt" $Database_outdir/"tmp.nc" $Database_outdir/"Regrid_fin.nc"
| grid_lcc.txt (1.17 KB) grid_lcc.txt | |||
| emissions.nc (7.42 KB) emissions.nc |
Replies (3)
RE: How to spatially distribute the values of a grid into another grid with different projection and resolution? - Added by Ralf Mueller about 5 years ago
hi Natalia!
I am not sure if CDO can fully understand your output grid description. You can test this by
cdo -f nc -topo,grid_lcc.txt topo_lcc.ncPlease check the output file wrt. the grid.
there is no indication of a land-sea-mask in your uploaded file. I found some 0-values in your field - if these are missing values, this should be explicitly said with defining the _FillValue attribute of your data variable to be 0. you can also do this with
cdo -setctomiss,0 emission.nc masked_emission.nc
hth
ralf
RE: How to spatially distribute the values of a grid into another grid with different projection and resolution? - Added by Natalia L about 5 years ago
Hi Ralf, thank you. Well, I didn't wrote the whole code where I used also a land mask. The command for the top does not work. It gets the error about the major axis , it finds it zero. This problem is also exists with the command griddes that's why I define the grit.txt by my own. The LCC has a custom spherical datum so a,b should be equal to 6370000.
I am writing below the whole code and I attach also the land use database I use with its grid description and the final output I get with remapycon.
#set missing value to 0
cdo setmisstoc,0 emissions.nc emissions_nozero.nc
#remap emissions files for the resolution of land use
cdo remapycon, "grid_reproj1.txt" emissions_nozero.nc "Regrid_lu_res_emissions.nc”
- Masking emissions based on lu
cdo ifthen land_use_map.nc "Regrid_lu_res_emissions.nc” "Regrid_lu_res_emissions_mask.nc”
#remap from the relosution of land use to the final grid lcc
cdo remapycon,"mygrid.txt" land_use_map.nc land_use_map_emiss1.nc
cdo setgridtype,curvilinear -setgrid,"mygrid.txt" land_use_map_emiss1.nc Regrid_fin1.nc
#set missing value to 0
cdo setmisstoc,0 Regrid_fin1.nc Regrid_output.nc
| grid_reproj1.txt (1.23 KB) grid_reproj1.txt | |||
| Regrid_output.nc (169 KB) Regrid_output.nc | |||
| land_use_map.nc (490 KB) land_use_map.nc |
RE: How to spatially distribute the values of a grid into another grid with different projection and resolution? - Added by Natalia L about 5 years ago
Natalia L wrote:
Hi Ralf, thank you. Well, I didn't wrote the whole code where I used also a land mask. The command for the top does not work. It gets the error about the major axis , it finds it zero. This problem is also exists with the command griddes that's why I define the grit.txt by my own. The LCC has a custom spherical datum so a,b should be equal to 6370000.
I am writing below the whole code and I attach also the land use database I use with its grid description and the final output I get with remapycon.
#set missing value to 0
cdo setmisstoc,0 emissions.nc emissions_nozero.nc#remap emissions files for the resolution of land use
cdo remapycon, "grid_reproj1.txt" emissions_nozero.nc "Regrid_lu_res_emissions.nc”
- Masking emissions based on lu
cdo ifthen land_use_map.nc "Regrid_lu_res_emissions.nc” "Regrid_lu_res_emissions_mask.nc”#remap from the relosution of land use to the final grid lcc
cdo remapycon,"grid_lcc.txt" land_use_map.nc land_use_map_emiss1.nc
cdo setgridtype,curvilinear -setgrid,"grid_lcc.txt" land_use_map_emiss1.nc Regrid_fin1.nc#set missing value to 0
cdo setmisstoc,0 Regrid_fin1.nc Regrid_output.nc